Analogue model of a FRW universe in Bose-Einstein condensates: 
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Analogue models of gravity have been motivated by the possibility of investigating phenomena not readily 
accessible in their cosmological counterparts. In this paper, we investigate the analogue of cosmological particle 
creation in a Friedmann-Robertson-Walker universe by numerically simulating a Bose-Einstein condensate 
with a time-dependent scattering length. In particular, we focus on a two-dimensional homogeneous condensate 
using the classical field method via the truncated Wigner approximation. We show that for various forms of the 
scaling function the particle production is consistent with the underlying theory in the long wavelength limit. In 
this context, we further discuss the implications of modified dispersion relations that arise from the microscopic 
theory of a weakly interacting Bose gas. 
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I. INTRODUCTION 

In the theory of quantum fields in classical backgrounds, 
some form of particle creation is typically expected when the 
metric is time-dependent; a commonly cited example of this 
is cosmological particle production in an expanding (or con- 
tracting) universe [1,2]. It is possible to simulate the analogy 
of this process in a Bose-Einstein condensate (BEC) when 
either: the external trapping frequency is time-dependent [3- 
7]; or the scattering length (within the low momentum ap- 
proximation of the two-body interaction potential) is time- 
dependent [8-10]. These treatments are often based on the 
acoustic {i.e., hydrodynamic) approximation where it is as- 
sumed that all excitations of the quantum field propagate as 
phonons with the same speed of sound. Moreover back- 
reaction of the excitations on the background field, and higher- 
order interactions between excitations, are neglected in the 
linearized theory. 

In general, the correct description of the dynamics of a BEC 
is a formidable problem due to the vastness of the Hilbert 
space, even for a system of just a few interacting atoms. In 
the lowest-order approximation, when all the bosonic atoms 
occupy a single quantum state, the ground state is well de- 
scribed by the Gross-Pitaevskii equation (GPE) — in this 
case, the field operator is replaced by a mean-field order pa- 
rameter Classical field methods (CFM) extend this formalism 
to include quantum fluctuations whereby the dynamics of a 
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multimode quantum field is approximated by the trajectories 
of classical variables in phase space. One such method is the 
truncated Wigner approximation (TWA), which is based on 
the Wigner representation of the density matrix. The TWA has 
been investigated by a number of authors [11-14] and more 
recently has been applied to a study of condensate collisions 
[15, 16]. 

In this paper we investigate the dynamics of a homoge- 
neous BEC in two spatial dimensions with a time-dependent 
scattering length, and in particular, map this problem to 
a Friedmann-Robertson-Walker (FRW) universe undergoing 
an expansion. We compare the analytically calculated parti- 
cle production from the acoustic approximation with the re- 
sults of numerical simulations based on the TWA. There are 
several benefits to this approach. Firstly, the TWA includes 
the effects of vacuum quantum fluctuations by sampling the 
Wigner distribution in the initial state. Secondly, the field dy- 
namics naturally include the nonlinear dispersion of the Bose 
system — in a cosmological context, this represents Lorentz 
symmetry breaking of the effective spacetime, which leads 
to necessary modifications of the standard hydrodynamic the- 
ory. Finally, the numerical simulations include the effects of 
back-reaction, which is difficult to otherwise include without 
resorting to higher-order methods such as the self-consistent 
Hartree-Fock-Bogoliubov approach. 

The outline of this paper is as follows: In Sec. II we show 
how the acoustic metric leads to an emergent FRW universe in 
a BEC, which leads to the prediction of quasiparticle produc- 
tion as discussed in Sec. III. In Sec. IV the these ideas are for- 
maUsed by the BogoUubov theory for a BEC. In Sees. V and 
VI respectively the preceding theory is used to quantitatively 
predict quasiparticle production in the acoustic approxima- 
tion, and in the more general case including "trans-phononic" 
modes, for a number of specific FRW universe models [17]. 
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Section VII introduces the TWA, which we use to simulate 
the dynamics of the BEC consistently with these scenarios. 
Section VIII presents the results of these simulations, which 
are discussed in Sec. IX. Finally, in Sec. X we conclude and 
discuss avenues for further work. 



n. EXPANDING UNIVERSE MODELS IN 
BOSE-EINSTEIN CONDENSATES 

Arguably, the most promising system for implementing 
analogue models of gravity is the BEC. This possibility was 
first considered by Garay et al. [18, 19] for acoustic black 
hole geometries, and further explored by Barcelo, Visser, and 
Liberati [20-22]. In particular, BECs have a number of de- 
sirable features that are necessary to be part of the analogue 
model programme, and more specifically for the subgroup of 
analogue models suitable for mimicking cosmological parti- 
cle production: 

(i) Hydrodynamics: In the long-wavelength limit, the 
mean-field equations of motion for a BEC take the form 
of classical hydrodynamics for a superfluid, which leads 
directly to the formulation of an emergent space time. 

(ii) Quantum theory: Bose-Einstein condensation in 
atomic vapours is a weakly interacting system, for 
which the microscopic quantum theory is well under- 
stood. The elementary excitations of the system are 
given by the Bogoliubov theory. The description for the 
excitations is valid beyond the hydrodynamic approx- 
imation, and incorporates "trans-phononic" physics, 
similar to many (not all) effective field theories, where 
"trans-Planckian" physics is beheved to break Lorentz 
invariance. 

(iii) Temperature: BECs in atomic vapours require tempera- 
tures close to absolute zero so that in principle it may 
be possible to observe cosmological particle produc- 
tion, without the presence of (larger) thermal fluctua- 
tions that would obscure the effect. 

(iv) Experimental advances: Recent experimental advances 
for the control of ultra-cold atoms mean that BECs can 
now be prepared and manipulated in many configura- 
tions. Notably, the use of magnetic and optical traps 
can lead to a variety of geometries, whereas using Fes- 
hbach resonances it is possible to vary the interaction 
strength between atoms, even by many orders of mag- 
nitude. With continued advances the experimental real- 
ization of analogue models should soon be achievable. 



A. Microscopic tlieory of the Bose gas 

Bose-Einstein condensation is characterised by a macro- 
scopic occupation of a single quantum state, typically the 
ground state of the system. This quantum degeneracy is 
achieved at very low temperatures where the de Broglie wave- 
length becomes comparable to the inter-particle spacing. In 



the second-quantised formalism, the effective Hamiltonian for 
a dilute Bose gas is 



H — Ho + Hi, 
where the single particle Hamiltonian is 



Hn 



dx (x) 



:^V2 + Kxt(x) 
Zm 



V'(x), 



and the interaction Hamiltonian is 



Hi^j I dx^t(x)V't(x)V'(x)V'(x). 



(1) 



(2) 



(3) 



Kxt(x) is any external potential (e.g., from a trap) and the two- 
body potential has been approximated via a contact potential 
U = Anfi^a/m in terms of the s-wave scattering length a, 
valid in the cold collision regime. The field operator ip{x.) 
annihilates a boson at position x and obeys the usual equal 
time commutation relations 

[V'(x,t), V;(x',t)] = [V>(x,0,7At(x',f)]=0, (4) 
[^(x,t), V^t(x',t)] = J(x-x'). 

The corresponding Heisenberg equation of motion for the 
field operator is 



ih 



av(x, t) 

'di 



2m 



^- (5) 



Without further approximations, this equation cannot be 
solved for a realistic system since the Hilbert space becomes 
prohibitively large, even for a system of just a few atoms. 

For the case of a weakly interacting Bose gas, and for 
low temperatures T ^ rcriticah a very useful approximation 
arises from the fact almost all the atoms reside in a single 
quantum state, so the system occupies only a fraction of the 
available quantum states. In this case it is useful to write 



V'(x, t) = Vlx, t) + (5(^(x, t), 



(6) 



where '(/'(x, t) = (V'(x, t)) is a mean-field term (the conden- 
sate wave function) and (5(^(x, t) is that part of the quantum 
field associated with quantum and thermal fluctuations, with 
((5(^(x, t)) =0. This is known as the Bogoliubov approxima- 
tion. This description itself may be treated with varying levels 
of approximation. In the very simplest approximation fluc- 
tuations are neglected altogether, and the resulting equation, 
known as the Gwss-Pitaevskii equation (GPE), is given by 



ih 



dt 



-^V2 + y,,,(x) + C/|^(x,t)|^ 



^(x,t).(7) 



This equation provides a description of the condensate in 
terms of a classical field. 
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B. Establishing the analogy 



It is now well established that effective spacetimes emerge 
from the microscopic theory of BECs in the long wavelength 
limit [20, 23]. In particular, the analogy between curved 
spacetime and BECs can be revealed starting from the mean- 
field description of condensates, given by the GPE (7). The 
complex order parameter t) can be written in terms of a 
real density and phase (i.e., the Madelung representation) as 



(8) 



Substituting this into (7) and equating real and imaginary parts 
leads to the pair of equations 



at m 



(9) 



and 



Linearizing about the fields n and 6 by setting 



n — * no + ni, 



equations (9) and (10) then lead to 

= -(711V6I0 +noV^i), 
ot in 



h—^ + uhi + —vOoVOi = 0, 

at m 



(11) 



(12) 



(13) 



where we have only retained terms up to first order in ni and 
9 1, and where we have defined the operator U by 



Uhi 



U - —D ' 
2m 



ni. 



(14) 



The differential operator D2 has been introduced, which to 
first order, takes the form: 



D2ni = -^?^o ^''^(V^V^)"-! + o ! — 



■(15) 



The operator D2 gives the first order correction from the in- 
clusion of the quantum pressure term. Identifying the back- 
ground velocity for the irrotational field as 



(16) 



and combining (12) and (13) yields a second-order differential 
equation for the perturbed phase, which can be written in the 



compact form [20, 23]: 

where we have introduced the symmetric matrix 



(17) 



with components [85] 



(18) 



r = 








= 








m 



(19) 



Note we have used the standard nomenclature where the 
Greek indices run from to d and the Roman indices from 
1 to d for d spatial dimensions. 

In the acoustic (or hydrodynamic) approximation the quan- 
tum pressure term is neglected. In this case D2 is set to zero 
(17) can be written more simply as 



1 







(20) 



where 5^;^ is the effective covariant metric tensor (with deter- 
minant g) given by 



.(x,i)=(-j 



(21) 



We have used the fact that the speed of sound in a condensate 
is given by 



= Uno/m. 



(22) 



The equation for the phase fluctuations (20) is formally anal- 
ogous to the dynamics of a massless and minimally coupled 
scalar field in a curved space-time [1]. It is worth emphasising 
that although the external potential V^^t does not explicitly ap- 
pear here, the field equation still depends on V^xt implicitly, as 
the the background geometry is determined by the stationary 
solutions of the GPE (7). 

In the acoustic approximation all collective excitations be- 
have as sound waves with the usual linear "relativistic" dis- 
persion relation 



LV = ck: 



(23) 



the quanta of excitations are thus phonons. An interesting con- 
sequence of Bogoliubov theory in Bose-Einstein condensates 
is that in general the excitation spectrum displays nonlinear 



4 



dispersion (see Sec. IV), being linear (i.e., phononic) for low 
|k| and becoming quadratic (i.e., free-particle like) at large 
|k|. When nonlinear dispersion is incorporated into analogue 
models of gravity it is equivalent to breaking Lorentz invari- 
ance [23-25]. We will return to this point in Sec. IIIB, and 
discuss the implications of this for the analogue model in a 
companion paper [17]. 

Equation (20) is the massless Klein-Gordon equation. The 
metric (21) has the signature ( — 1-++) and so the effective 
spacetime represents a Lorentzian geometry. It is interesting 
to note that an effective "relativistic" wave equation (for exci- 
tations in the condensate) arises from the equations of motion 
for a non-relativistic fluid [26, 27]. 

The hydrodynamic (i.e., long wavelength) description of a 
BEC is commonly believed to lead to an effective field the- 
ory (EFT) for curved spacetime in the same sense that clas- 
sical gravity is expected to emerge as the low-energy EFT of 
a full theory of quantum gravity [28]. (The analogue grav- 
ity EFT arising from a BEC does not impose a strict cut-off 
and so includes so-called trans -phononic physics [25, 29-32] 
— we return to this point in Sec. IIIB). However, it should 
be noted that the emergent spacetime is itself embedded in an 
low-energy EFT for the Bose system; here the cut-off is de- 
termined by the validity of replacing the two-body interaction 
potential by an effective potential that is characterised by the 
s-wave scattering length. In practice, the cut-off is enforced 
by the the use of a projector in mode space. 

We further emphasise that the above analogy only holds for 
massless spin-0 particles; in general, it is possible to mod- 
ify the formalism to include massive modes at the expense 
of dealing with a more complex configuration, e.g., a two- 
component BEC [25, 29-34]. Moreover, it is possible to de- 
velop an analogue model of gravity for spin-1 fields, e.g., in 
dielectric media [35]. 



C. FRW analogue model 

Time dependence can enter the metric tensor (21) in any of 
the parameters no, c, and v. We focus on the case where the 
background flow is zero (v = 0), and the system is homoge- 
neous with a density no that is constant through space. With 
this choice of parameters the effective metric (21) becomes 



'no\ — 



c2 ; 



: 6,. 



(24) 



This represents a Friedmann-Robertson-Walker spacetime 
[8, 9]. (Technically, a fc = spatially flat FRW space- 
time.) Such spacetime geometries are conformally flat, and 
at any particular time the spatial geometry is simply that of 
Euclidean 3-space. The time dependence is contained entirely 
in the speed of sound given by 



c{tf 



U{t)na 47rft2 



m 



noa(t), 



(25) 



with atoms of mass m, scattering length a and number den- 
sity n. We introduce the dimensionless scaling function b{t) 
so that the interaction strength (or equivalently the scattering 
length) has the time dependence 



U{t) = Uob{t), 



(26) 



where Uq = U{to) at an initial time to; therefore we take 
b{to) = 1. From (25) the time dependence of the speed of 
sound is thus 



c(t) = Co b{t) 



1/2 



(27) 



In practice a variation in the interaction strength is possible 
by using a Feshbach resonance [8, 36, 37]. If b{t) is decreas- 
ing with time we have an expanding universe model, whereas 
if b{t) increases with time we have a contracting universe 
model. It is useful to introduce X the expansion of the uni- 
verse between two times to and tf in the following way: 



X = bitf)"- 



(28) 



Please note that one can always without loss of generality 
choose 6(to) = 1. The line element for the FRW universe 
we have described is given by 

dslff = nl [^clb{trdt^ + 6(t)"-idx2] , (29) 

where we have introduced the conformal factor 



(30) 



f^o("-o,co,cf) = — 



Co 



which is independent of space and time, as well as the 
dimension-dependent exponent 



d-2 
a ~ — . 

d- 1 

Hence, in c? = 2 spatial dimensions we get 

X = b{tf)-\ 



(31) 



(32) 



as the overall expansion of the universe between Iq and t f . 

In what respect do we have an expanding universe, given 
that the condensate is contained in a physically fixed volume 
V7 A decrease in the scattering length corresponds to a de- 
crease in the the speed of sound propagating in the conden- 
sate; therefore any acoustic excitations will propagate with 
decreasing speed in the condensate as time passes. To an ob- 
server at rest in the effective spacetime, a decrease of the speed 
of sound is thus indistinguishable to an isotropic expansion of 
the spatial dimensions. In this sense, (24) corresponds to the 
notion of a spacetime undergoing an effective expansion. 

It is not straightforward to define either an apparent or event 
horizon in the model considered here as the system is homo- 
geneous, and the background velocity is therefore the same 
(i.e., zero) everywhere. This is further complicated by the fact 
that the causal structure of the effective spacetime should be 
determined by the maximum signal velocity (i.e., group ve- 
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locity), which is effectively infinite here owing to the super- 
phononic modes in a BEC. Alternatively, analogue models 
where the background velocity depends on the radial position 
[e.g., when the trapping potential is switched off and the con- 
densate is free to expand) have been studied in [3-7]. In the 
present situation, where the scale factor aFRw(0 or equiva- 
lently b(t) contains all the geometric structure for the space- 
time, it is necessary to perform the usual analysis in terms of 
cosmological horizons to determine the overall causal struc- 
ture of the spacetime [38]. 



D. Two dimensional model 

To facilitate the numerical calculations required by the clas- 
sical field simulations that we present in Sec. VIII, we con- 
tinue with d ~2 spatial dimensions. The reduced mode space 
for d ~ 2 greatly decreases the computation time, but still 
leads to a satisfactory description of the system. In particular, 
we expect the extension of the numerical simulations to d = 3 
to lead to qualitatively similar results [86]. Moreover, while it 
has been shown [6, 7] that for d = 2,a condensate undergoing 
free expansion leads to a scalar field equation that does not de- 
pend on the scaling factor in co-moving coordinates [87] — so 
that in this case there can be no particle production — as we 
will see, the present model for a 2 + 1 dimensional FRW uni- 
verse leads to a scalar field equation that does allow particle 
production. 

For 2 + 1 spacetime dimensions, a{2) = and the line- 
element simphfies to 



]dt^ + b{t)-^dx^] . (33) 



A further time transformation is not required as laboratory 
time and proper time (for a comoving observer) are of the 
same form. Although (33) includes a time- and space- 
independent conformal factor, this does not affect the dynam- 
ics of the field (in this case the mode functions need to be 
normalized for consistency with the Bogoliubov theory). The 
scaling factor aFRw(i) for a FRW universe that is familiar from 
cosmology is related to b{t) by 



aFRw(i) = b{t)'^^^ 



(34) 



for d = 2. (We always explicitly specify the FRW subscript 
for the cosmological scale factor so that this quantity is not 
confused with the s-wave scattering length a.) 

The reduction of the model to two dimensions requires a 
modification to the nonlinearity that appears in the GPE (7), 
and therefore also the resulting field equation (17). To see this, 
we assume the transverse z dimension is tightly confined with 
the trap lengths satisfying L, ^ L^, Ly, and further that the 
~ ^ for the transverse dimension where ^ is the healing 
length of the condensate. The scattering is still determined by 
the three-dimensional scattering length so that this is called 
a quasi-two dimensional geometry. However, the system re- 
mains in the ground state of the transverse dimension because 



the energy required for transverse excitations is much larger 
than for longitudinal excitations. The wavefunction is then 
separable as \E'(x, i) = ijj{x,y,t)C,{z). Assuming the con- 
densate is homogeneous in the z direction, the GPE can be 
rewritten as 

dip{x,y,t) ( ^ \ 

dt " \2^^ +V,nj4,{x,y,t) 

+U2v\^{x,y,t)\''4,{x,y,t), (35) 

where the effective nonlinearity is U2D — U/Lz- This does 
not affect the form of the resulting calculations for particle 
production, but should be noted when we determine suitable 
parameters for our simulations in Sec. VII F. (Note that be- 
cause we are changing the atomic interactions by several or- 
ders of magnitude during the proposed experiment, this will 
also affect the transverse trapping.) 



III. FIELD QUANTIZATION AND PARTICLE 
PRODUCTION 

The field 6{t) is quantized using the plane wave mode ex- 
pansion: 



(36) 



where &k and b\^ are the annihilation and creation operators 
respectively for the quasiparticle modes. 

In flat (Minkowski) spacetime we can associate the positive 
and negative frequency solutions of (20) with annihilation and 
creation operators respectively. In curved spacetime, this as- 
sociation is not always possible as different observers experi- 
ence different vacua; in the language of general relativity, the 
spacetime need not have a (globally defined) time-like Killing 
vector field so that the positive frequency solution is not nec- 
essarily unique and typically is not an eigenfunction of dt- A 
consequence of this is that the "natural" choice of a Fock vac- 
uum according to 5|0) =0 depends in general on the choice 
of coordinates; that is, the measurement of particle content in 
curved spacetime is said to be "observer dependent". 

The calculation of particle production follows the standard 
methodology [1, 2]. We define in and out regions respectively 
as asymptotically flat regions with t — > — cx) and t +00. 
(While the existence of asymptotically flat regions cannot be 
assumed for cosmological models, it is certainly possible to 
emulate this scenario in BEC experiments.) We can write a 
mode expansion (36) for the in and out regions in terms of 
mode functions or X\t^ respectively. Because both sets 
of modes are a basis set for the field, they are related by the 
Bogoliubov transformation 



, .out m I ID ,,in* 

Xk - "kXk +P-kX-k- 



(37) 



Note the spacetime has translational invariance as a symme- 
try so the field (36) can be expanded in the same set of plane 
wave modes e''' for both the in and out regions, and there- 
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fore /3_k = /3k also. By convention the mode functions are 
normalised according to the Klein-Gordon inner product. In 
momentum space this takes the form 

-«M^[Xk,Xk] = l, (38) 

where W[fi,f2] = fi{dtf2) - {dtfi)f2 is the Wronskian of 
two functions fi and /2. 

For each asymptotic region the vacuum state is defined by 
the requirement 

6if |0-) = 0, 

^OUt|QOUt^ = 0. (39) 

The in and out mode operators are also related by 

bw' = "k^L" + /3-kfotJk- (40) 
The number of particles created is then given by the quantity 

iV^* = (0'"|fe[°^*C*|0'") 

= |/?k|'. (41) 

The standard procedure to calculate /3k involves solving the 
field equation (17) for the in and out regions by using (37) 
with the appropriate boundary conditions. 

It should be stressed that here particle production refers to 
the production of (massless) quasiparticle excitations in the 
Bogoliubov basis, that approximately diagonalizes the many 
body Hamiltonian ( 1 ) for the Rose gas. As (40) couples modes 
of momenta k and — k this process is often referred to as pair 
production and is associated with the formation of squeezed 
states [39-43]. 



1. Choosing a Fock vacuum 

When the in and out regions of the expansion are not 
asymptotically flat, a vacuum state cannot be unambiguously 
defined for either case. However, clearly the procedure to cal- 
culate particle production outlined above requires a choice of 
Fock vacuum, and moreover, that choice should lead to phys- 
ically reasonable results. We therefore mention two choices 
that have been developed to deal with this situation (although 
there are many more): 

The instantaneous Minkowski vacuum: this is the 
state that corresponds to the instantaneous diagonaliza- 
tion of the Hamiltonian at a given time (and is there- 
fore also the state that minimises the energy). (This 
procedure is also known as Hamiltonian diagonaliza- 
tion.) This choice may be problematic in some cases 
as it can lead to situations of infinite particle production 
even though the expansion may be smooth and finite 
[2]. 

The adiabatic vacuum: this approximation can be 
used for modes that experience a sufficiently slow ex- 
pansion and is formulated in terms of the WKB approx- 



imation [1, 44]. However, the adiabaticity requirement 
means it has limited applicability. 

In Sec. V C we calculate particle production for the case of a 
de Sitter expansion, for which there are no static regions — 
in this case we therefore utilize the instantaneous Minkowski 
vacuum to calculate the Bogoliubov coefficients at a given 
time during the expansion. Once the expansion has stopped 
the final Bogoliubov coefficients unambiguously correspond 
to physical particle production. The problem of infinite parti- 
cle production does not occur in this case, and as we shall see 
from the classical field simulation results in Sec. VIII, particle 
production is further suppressed for short wavelength modes 
because of super-phononic dispersion. 

2. Commutation relations 

In addition to identifying an effective metric required to 
make the analogy with quantum field theory in curved space- 
time, another prerequisite is that in quantizing the field the 
operators should be annihilation and creation operators in the 
usual sense. That is, they should obey the correct commuta- 
tion relations so that we can define quanta of the field in the 
Fock basis. This has indeed shown to be the case for the lin- 
ear excitations of a BBC [45, 46], and also follows from the 
discussion in Sec. IV A. 



3. Relevant scales 

There are two relevant scales for particle production in a 
FRW-type analogue model: 

1. Hubble parameter - Following [8] for the metric (33) 
the Hubble parameter is given, for d = 2, in terms of 
the scaling function h{t) in laboratory time by 

H^^^ = -\\. (42) 
apRw ^ 

This quantity, which is time-dependent in general, cor- 
responds to the rate of expansion for the universe. If 
H > LUk for a mode frequency labeled by k the dynam- 
ics are Hubble-dominated and we expect the mode to 
be to be non-oscillating, whereas H ^ lo^. implies the 
mode is oscillating and in the adiabatic regime. In the 
latter case this is what cosmologists refer to as a para- 
metrically excited mode. Clearly a large value for H 
is favourable to particle production — however H can 
not be made arbitrarily large as the approximations that 
lead to the effective field theory would then be violated 
[8, 47]. 

2. Healing length - In a BEC, the healing length is a dis- 
tance over which localized perturbations in the conden- 
sate tend to smooth out, and is given by 

^^^^iobitr'/'. (43) 
V Zmc 
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That is, if we define the cross-over from phonon to free- 
particle behaviour for the BogoHubov spectrum (62) as 
h^k1/2m = Un then kc = 1/^. Modes for which 
fc ^ 1/^ correspond to collective excitations of the 
condensate (phonons) and couple to the effective time- 
dependent curved spacetime whereas modes with k ^ 
1 are particle-like and are relatively unaffected by the 
emergent spacetime geometry. Alternatively stated, the 
spacetime appears locally flat and time-independent to 
modes with sufficiently short wavelengths. 

Particle production into a mode k then proceeds as fol- 
lows: During expansion, while H > ujk and k < 1/^, the 
mode evolves non-trivially and particle production will cer- 
tainly occur As the expansion proceeds, particle production 
can "switch off" for two reasons: Either the healing length in- 
creases until the mode becomes particle-like and particle pro- 
duction slows, ceasing altogether when fc 1/^; or alterna- 
tively, if the expansion slows, then after some time H < Uk 
so that the field begins to oscillate relatively freely, and addi- 
tional particle production also ceases. In either case, the mode 
no longer evolves, and the occupation number of the mode be- 
comes constant. 

In general, it can be difficult or impossible to analytically 
solve the field equation (17) for the mode functions. More- 
over, if the in I out region is not flat, there is no preferred 
choice for the initial / final vacuum state. However, when the 
field evolution is sufficiently slow — i.e., adiabatic — the 
WKB approximation can be used to calculate particle pro- 
duction (see [44] and references therein). It is worth noting 
that this approach leads to a Planckian spectrum in the lowest 
order approximation [1]. 



4. A note on freezing 

For our analogue model of an expanding FRW universe, in 
the acoustic approximation we have LOk{t) ~ c{t)k, which 
decreases as the expansion proceeds (this is equivalent to the 
usual notion of cosmological redshifting of modes that occurs 
during inflation). This means that depending on the form of 
the expansion, a mode that is initially oscillating {H ujk) 
may enter a Hubble-dominated era {H ^ ujk) after a suffi- 
ciently long expansion — this is the mechanism for freezing 
of modes that is familiar from inflationary cosmology. The 
tenn freezing does not necessarily imply that particle produc- 
tion ceases, but rather, that the mode no longer oscillates. 

The situation is quite different when nonlinear super- 
phononic dispersion is included — in this case, the healing 
length ^ also decreases for increasing time, which means that 
a mode that was initially phononic will crossover into the adi- 
abatic free-particle regime for a sufficiently long expansion. 
In this case therefore, the notion of freezing can only occur 
in a transient regime that depends on the mode wave-vector k 
and the form of the scaling function b{t) (and therefore on the 
Hubble parameter also). A more detailed discussion can be 
found in [17]. 



A. Acoustic approximation 

In the acoustic approximation, the correspondence between 
the massless minimally coupled scalar field equation for a 
FRW universe and the equations of motion for linearized fluc- 
tuations in the EEC is exact. The dynamics of the low momen- 
tum modes are expected to fall into this regime, and therefore 
it is of interest to examine this case first. 



1. Conformal time 

We define a coordinate transformation from laboratory time 
t to conformal time ?/ by 

dri = ^/b(J)dt. (44) 
The line element (33) then reads 

dsls = nl b{r])-' [-cldrf + dx2] , (45) 

In (2 + 1) dimensions the equation of motion for the field 
obtained from (17) in conformal time becomes 

^'^'-lTWf^'-'"^^"' = ''- ^^^^ 

Note that the coefficient of djj9 is the Hubble parameter given 
by (42). The calculation of particle production requires the so- 
lution of this field equation for the mode functions xiv)- This 
task is assisted by reducing the field equation to standard form 
where the first order derivative of the field operator does not 
appear; this can be achieved (for instance) by introducing an 
auxiliary field of the form ^{rj) = b{rj)0, showing mathemati- 
cally equivalent dynamics to 0. Alternatively we can consider 
the transformation to auxiliary time as follows. 



2. Auxiliary time 

Defining the auxiliary factor by 

A = ny^b[t)-^l'^, (47) 
and the auxiliary time by 

the line element (33) transforms as 

dlls = -A^dP + A^dx^ . (49) 

For a massless scalar field, the corresponding field equation in 
(2 + 1) dimensions is 

dfe - A^V^e = 0. (50) 
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Using the field mode expansion given by (36) we then get a 
time-dependent harmonic oscillator for each mode 

9?Xk + c:'kWxk = o (51) 

where <D^(t) = K^k"^ = ilo^(i)^^fc^ is the oscillator fre- 
quency. While the auxiliary time approach is useful in the 
acoustic approximation, it does not generally lead to the sim- 
ple form (51) when the quantum pressure term is included (see 
Sec. IIIB). 

B. Beyond the acoustic approximation 

The description from the previous section was given within 
the acoustic approximation, valid only for long wavelength 
modes, whereby the quantum pressure term is omitted from 
the field equation (17) and the resulting equations for the FRW 
analogue model given by (46) or (51). We can extend the 
analysis to higher momentum modes by including the quan- 
tum pressure term; this leads to a super-phononic type dis- 
persion relation for high frequency modes and an appropri- 
ately modified field equation as we presently discuss. More- 
over, this description corresponds to the presence of "trans- 
Planckian" like modes, and leads to the idea of Lorentz viola- 
tion. Such modified dispersion relations have similarly been 
incorporated into some studies of inflationary cosmology — 
we discuss this here briefly, but dedicate our companion pa- 
per [17] to exploring this subject. In Sec. VI we will consider 
two specific forms of the scaling function in this regime: the 
limiting case of a sudden transition (Sec. VIA) and a cyclic 
universe model (Sec. VI B), both which include the effect of 
quantum pressure. 

1. Nonlinear dispersion 

We include the quantum pressure term by using (77) so that 
the field equation (17) in momentum space becomes 

dtXk - \^j^dtb{t) dtXk + 0Jk{tf Xk = 0, (52) 

where we have now defined 

^.(i)' = |^(^^ + 2C/Wnoj. (53) 

It should come as no surprise that we have recovered the Bo- 
goliubov dispersion relation for a weakly interacting Bose 
gas (see Sec. IV below). In particular, we see that if we set 
h{t) ~ const, then the field equation describes the dynamics of 
each Bogoliubov mode for a time-independent Hamiltonian. 

In general, it is difficult to solve (52) for the mode func- 
tions, and to hence calculate particle production. This situa- 
tion often persists even though it may be possible to reduce 
(52) to standard form where the first order derivative term 
dtOk does not appear. We can, however, make a qualitative 
statement about particle production for high momenta modes 



by considering the crossover from phononic to free-particle 
modes as determined by the healing length (43). For modes 
that satisfy k 3> fee, the field equation takes the approximate 
form 

where 



represents the single particle energy for a non-interacting gas. 
That is, the time dependence from b{t) above is largely sup- 
pressed, and each mode with k ^ kc evolves trivially as a 
time-independent harmonic oscillator, remaining in its initial 
vacuum state. Therefore these modes experience no particle 
production. 



2. Lorentz violation 

Planck-scale Lorentz violation is a feature of some theo- 
ries of quantum gravity that can be modeled by the presence 
of a modified dispersion relation at the Planck scale [25, 29- 
32, 48]. In this sense, for a BEC the healing length ^ provides 
the analogue of the Planck scale, characterising the cross-over 
from phononic (fc <C 1/0 to free-particle (fc ^ like 
modes. Specifically, the field equation (52) is not Lorentz in- 
variant because of the nonlinear super-phononic dispersion of 
the modes (53) at large momenta (this clearly evident from 
Eq. (54) where the second term varies as ~ k'^). On the other 
hand, for small momenta {i.e., in the acoustic approximation) 
the dispersion relation is linear in k and the field equation re- 
duces to the Lorentz invariant form (46). It should be stressed 
however, that the effective quantum field theory for the BEC is 
still valid as long as atomic interactions can be characterised 
by the s-wave scattering length, which is true in general for 
some cut-off in wave-vector fccut-off 

We note that modified dispersion relations have sometimes 
been used to study the "trans-Planckian" problem in inflation- 
ary cosmology [49-52] — the results there showed that cer- 
tain modifications to the dispersion relation could lead to sig- 
nificant deviations for the density fluctuation spectrum when 
compared to the unmodified dispersion relation {i.e., the usual 
model of inflation). 



IV. CONNECTION WITH BOGOLIUBOV THEORY 

Thus far, we have derived a field equation for a scalar field 
propagating in an effective spacetime. The quanta of this 
scalar field must correspond to the linearised quantised ex- 
citations {i.e., phonons) of the quantum field for the theory 
to be consistent. Therefore, at this point, to make this con- 
nection explicit, it is worth reviewing the theory of quantum 
excitations in BECs, which is well described by Bogoliubov's 
theory of excitations for a weakly interacting system. 
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A. Bogoliubov Theory 



where 



(64) 



While the GPE has been extremely successful for describ- 
ing mean-field effects {i.e., classical dynamics), it neglects 
quantum and thermal fluctuations. To first order, we can in- 
clude quantum fluctuations in this description by considering 
the theory of elementary excitations for a weakly interacting 
Bose gas, first formulated by Bogoliubov for the homoge- 
neous case [53]. The standard procedure is to first expand 
the field operator in a plane-wave basis 



H^, t) 



,'ikx 



ak(i)- 



(56) 



If the number of atoms in the condensate is large {i.e., the 
field is highly condensed) then replacing ao, a\ \fN^ and 
retaining terms of at least order A^o. the Hamiltonian (1) can 
be approximately diagonahsed using the Bogohubov transfor- 
mation 



(57) 



when Wk and Vk are chosen so that the new operators 6k and 
satisfy the commutation relations for Bose field operators 



2 2 1 
"k - «k = 1- 

This leads to the result that [54, 55] 

1 Ak 



Mk 



Wk 



(58) 



(59) 



with 



A. 



1 

Unn 



-{(i + Uon) + ^el{el + 2Um) , (60) 
where n = Nq/V is, the density. The resulting Hamiltonian is 
H « Hfiog = ^0 + ^ (J^lh^, (61) 

k^^O 

with a constant Eq and where the quasiparticle excitations 
have the energy spectrum 



Ck 



(62) 



in terms of the single particle energy for a non-interacting 
gas, see Eq. (55). In the Bogoliubov approximation, explicitly 
pulling out a phase factor depending on the chemical poten- 
tial ji = J/riQ. the Bose field operator can conveniently be 
expanded as 



V'(x,i) = e-'''*/''[V^+5(^(x,t)] 



(63) 



Note that the time-dependence of each mode is fully con- 
tained in the mode functions Uk{x, t) and Va;(x, t). When the 
Hamiltonian is time-independent {i.e., U constant) the time- 
dependence is purely oscillatory and the number of quasipar- 
ticles in each mode (blbk) is a constant of the motion. 

Alternatively when we linearize the density-phase represen- 
tation of the field operator we find 



If hi and 6i are Hermitian we can write 



1 



no I 



(Sip — Sip^). 



(65) 



(66) 



(67) 



The commutation relations for the field operator (4) can be 
used to show that the density and phase fluctuations satisfy 
the commutation relation 



[ri,i(x),^i(x')] = i6{-K - x'). 



(68) 



Using the mode expansion (36) we expand the field 6i in 
Fourier (plane-wave) modes as 

^1 - 7^ E[^^'''"Xk6k(0 + e-^'^-xlblit)], (69) 



and similarly for hi 



ni = ^ 2^[e**'-"nk6k(i) + e^'^^-^^^t (t)]- (70) 

^ k 



The commutation relation (68) thus reduces to: 



n-kXl - nlXk = i- (71) 



We can expand S(p, hi and di in the same set of plane wave 
modes, using the mode expansions (69) and (70) and 



and 

so that the Fourier components are 

Uk{t) = ^ — nk{t) + iy^Xk{t), 



(72) 



(73) 



(74) 
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Vkit) 



nk{t) - iy/ri^Xk{t)- 



(-75-) The mode functions are then given by 



Clearly these mode functions are consistent with the require- 
ment that the Bogoliubov modes are normalised by jiifcP — 
= 1. We note that these are general expressions, which 
are vaUd within the linearized theory of excitations regardless 
of what form the mode functions take. That is, they are vahd 
for arbitrary forms of the scaling function b{t) so long as the 
mode functions nk{t) and Xk{t) can be found. 



Xk 



U 
2hujk 



2U 



(82) 



(83) 



Using the general expressions (74) and (75) the mode func- 
tions are given by 



B. Bogoliubov modes — Minkowski spacetime 

For the homogeneous model presented here the differential 
term (15), accounting for the quantum pressure term, takes the 
simple form 



1 2- 
D2ni = V n\. 

2rio 



(76) 



Rearranging (13) and extracting the Fourier component then 

gives 



hujk 


j2Uno 




V f^k 


huk 


hUno 


2Uno 


V fi^k 



(84) 
(85) 



h 



2Unr 



dtX-k- 



(77) 



D. Quasiparticle production 



Using (77) with the mode expansions (69) and (70), the com- 
mutation relation (68) is satisfied when the mode functions 
take the form 




(78) 



These are the positive-frequency solutions for a time- 
independent Hamiltonian — that is, when when U is constant 
— and correspond to an effective spacetime geometry that is 
Minkowski flat. The mode functions (74) and (75) can there- 
fore be written 



Uk 



(79) 
(80) 



The concept of particle production in our analogue model 
can be made explicit in terms of the Bogoliubov theory out- 
lined above. To calculate the quasiparticle number in each 
mode after a finite expansion of duration t it is necessary to 
project onto the Bogoliubov basis that instantaneously diag- 
onalizes the many-body Bose Hamiltonian at t. (That is, we 
take the instantaneous Minkowski vacuum as the zero parti- 
cle state.) In what follows we take the initial condition as the 
quasiparticle vacuum at t = so that 6k(0)|0) = and there- 
fore Nk{0) — 0. Here the mode functions are determined by 
the solutions to the time-independent case with U (t) = Uq, 
which yields the mode functions for a Minkowski spacetime. 
With this premise, we now proceed to calculate the particle 
production for each mode for a given expansion. The Bo- 
goliubov theory of a weakly interacting Bose gas predicts a 
non-zero depletion even at zero temperature; the real particle 
annihilation operator is given by the time-dependent canonical 
transformation 



ak(t)=t*r(t)6k(0)+t'r*(t)6Lk(0) 



(86) 



C. Acoustic modes — Minkowski spacetime 

To facilitate the computation of particle production we 
again consider the acoustic approximation, which is valid for 
low momenta when fv^k"^ /2m ^ i/ng so that the quantum 
pressure term can be neglected. In this case (77) reduces to 



nk 



U 



dtX-k 



(81) 



where uj^^ and v^!^^ are solutions for the mode functions dur- 
ing the expansion — these must coincide with Minkowski 
mode functions at t = with U{t) = Uq. The projection 
into the Bogoliubov basis at t is 



(87) 



where u™* and u^"* 



are given using (74) and (75) and the 
Minkowski mode solutions (82) and (83) with U{t) = Uq/X 
in terms of the expansion X. The particle production in each 
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mode at time t is then given by 

, out* / 



\ur*it)vr*{t)-vr*it)ur*it)f 
1 

-2\urmvnt)\ 



\nnt)\'+no\xnt)\' 



1 

4no 



\nnt)\'-no\xr{t)\' 



1 



:-(88) 



where we have used (74) and (75). Evidently, the central task 
is to solve the field equation (20) for the mode functions u™^ 
and v'^^ for a given expansion b{t). This procedure will be 
applied in Sec. V C to calculate the quasiparticle production 
for the case of de Sitter expansion (in the acoustic approxima- 
tion), and in Sec. VI A for the case of sudden expansion (with 
high-frequency dispersion). 



by applying the boundary conditions from (37) and its first 
derivative for the Minkowski in and out modes at to = 0. We 
find 



and 



so that 



l/?k 



Further noting for an expansion X that (I'£"'/'^k' 
the particle production is then 



(91) 



(92) 



(93) 



i/Vx, 



(94) 



V. QUASIPARTICLE PRODUCTION IN ACOUSTIC 
APPROXIMATION 

We presently provide analytic solutions for quasiparticle 
production in the acoustic approximation for three different 
expansion scenarios: (i) Sudden transition; (ii) tanh expan- 
sion; and (iii) de Sitter expansion. Scenarios (i) and (ii) both 
have asymptoticly flat in and out regions, so the calculation of 
particle production follows the standard procedure from sec- 
tion 111. Scenario (iii) however does not have asymptoticly flat 
in and out regions for a finite time expansion, and therefore we 
must resort to the Bogoliubov theory of the previous section 
to calculate the particle production. 



A. Sudden transition (2 + 1 dimensions) 

A simple example of particle production is given by the lim- 
iting case of a sudden expansion. Here, the interaction is in- 
stantaneously switched from U toU/X at some time which 
we take as fo = for convenience. The scaling function is 
given by 



This yields, for example: 7V£"* « 0.4 for X = 10; iV£"* « 2 
for X = 100; and 7V^"* « 11 for X = 2000. 

Equation (94) provides an upper limit for the particle pro- 
duction in each mode [56]. This quantity does not depend on 
mode number — a feature that reflects the fact that all modes 
experience a sudden change in the effective spacetime geom- 
etry. We shall see in section 111 B that including the quantum 
pressure term (i.e., nonlinear dispersion) in our formulation 
leads to suppressed particle production for increasing |k|. 

It should be noted that the sudden transition corresponds to 
a delta function for the Hubble parameter iJ at t = 0; this is 
physically unfeasible as any change in the s-wave scattering 
length via a Feshbach resonance would require a finite time 
in practice; moreover, a very rapid change in the scattering 
length is not possible since the low momentum approximation 
for the T-matrix scattering potential is no longer valid [8, 47]. 
In spite of this, the sudden transition still provides a useful 
prediction for comparison with the results of the classical field 
simulations. 



B. tanh expansion (2 + 1 dimensions) 



b{i) 



1 



(89) 



where H is the Heaviside step function. The case of parti- 
cle production for a sudden transition has been previously ex- 
plored by Jacobson for a parametric oscillator [41]. 

Using the normalisation condition (38) the positive fre- 
quency solutions to (51) are given by 



./ out 



1 



. . ~ in/ out 

2^k 



(90) 



It is straightforward to calculate the Bogoliubov coefficients 



One non-trivial form of the metric tensor for which the par- 
ticle production can be calculated analytically is the case of a 
tanh function expansion, with asymptotically flat in and out 
regions. This was first considered by Bernard and Duncan 
[57], and Birrell and Davies [1] for a massive scalar field, and 
then by Barcelo et al. [8] for a massless scalar field. Similarly 
to [8] we consider the case of tanh expansion, but for 2h-1 
dimensions. 

In particular (51) can be solved exactly when the auxiliary 
factor A has the time-dependence 



A? 



A3 



A3 



tanh ( 

2 \U 



(95) 
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for some time constant tg that determines the rate of expan- 
sion. Noting that 



1 



(96) 



the scaling function with respect to auxiliary time is 

'A3 + A3 A3„A3 



bit) 



= 2 



2 2 

1 + X + {X -1) tanh 



tanh 



(97) 



for an expansion X; we have implicitly assumed b{ti) = 1. 
We can also write 

Co J bit) 



2col 



il+X)i+ iX-l)i,log 



cosh 



(98) 



This is not easy to invert but b{i) and t{i) define a parametric 
curve for b{t) — this relation is required for implementing 
a tanh expansion in laboratory time simulations. With the 
conformal factor given by (95) it is possible to calculate the 
particle production in each mode exactly. Using the result in 
[8] we get 



^out 



sinh^ 


TrkisQaiV 


X-l)/2 




sinh [Trfcisfio] sinh 




X 



(99) 



Note, in the limit is ^ this reduces to the sudden transition 
result (93) as expected. 



C. de Sitter universe (2 + 1 dimensions) 



The case of the de Sitter universe is particularly relevant 
in cosmology. The inflationary model of the early universe is 
thought to include a de Sitter phase of rapid expansion which 
ultimately accounts for the inhomogeneities observed in the 
present universe [58]. The de Sitter spacetime is a solution 
to the Einstein's field equations with a positive cosmologi- 
cal constant, and has a high degree of symmetry. It has been 
shown that an observer moving in a time-like geodesic will 
measure a thermal spectrum — this particle production from 
cosmological horizons being related to the Hawking and Un- 
ruh effect (for event and particle horizons respectively). This 
result was first derived by Gibbons and Hawking using the 
path integral formalism [59] and has been subsequently veri- 
fied by applying the method of Bogoliubov mode mixing [60] . 



m 



Region I 



Region II 



Region III 



to 



FIG. 1: Schematic of de Sitter expansion for the FRW analogue 
model; for region I (t < to), there is no expansion and the mode so- 
lutions are those of a Minkowski spacetime with U = Uo', for region 
II (to < t < tf), there is a de Sitter type expansion and the mode 
solutions are non-trivial; finally for region III (t > tf) the expan- 
sion is turned off and the mode solutions are those of a Minkowski 
spacetime with (7 — Uo/X. 



1. Scaling function 

To map the FRW analogue model to a de Sitter spacetime, 
the scaling function for the scattering length in laboratory time 
(which is equivalent to proper time for two dimensions) is of 
the form 



bit) 



-t/ts 



(100) 



with the scaling unit ts that determines the rate of expansion. 
In this case the Hubble parameter (42) is given by i7 = 1 /2ts. 

We further consider the transformation to so-called confor- 
mal time (denoted rj) by drj ~ \/bit) dt. In this case we get 



ry = -2t,e-*/2*= 



t > 0. 



(101) 



The following limits are evident: (i) rj = —2ts for t — 0; and 
(ii) 77 ^ as < ^ +00. We also have: 



biv) = 



-2ts <r]<0. 



(102) 



2. Mode solutions 

The field equation (46) then yields the second-order differ- 
ential equation for the mode functions 



d^Xk 19xk , 2,2 n 
-jr^ ^ ^ Cq/c Xk = 0. 



(103) 



This is a Bessel equation and the solution [61, Eq. 9.1.52] is 
given in terms of Bessel functions of the first and second kind 
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[88] 

Xk^AkV Ji{-LJk{0)v) + Bk V Yi{~uJk{0)v) (104) 

for some undetermined constants Ak and Bk', we have defined 
the frequency u!k{Q) ~ cok at t = 0. Using (81) and [61, 
Eq. 9.1.27] it can be shown that the density fluctuation mode 
function is given by 



hLUk{0) 

Uk = -2ts — — — [AkJa{-uJkil) 



BkY^{~ujkTl)].{lQ5) 



With no loss of generality we choose h{t) = 1 for t < — 
i.e., for 77 < -~2ts. Thus i = (or 77 = — 2ts) corresponds 
to the end of the in region, whereas t > tf corresponds to the 
out region (see Fig. 1). To determine the particle production it 
will be necessary to calculate the Bogoliubov coefficients ak 
and (3k between these two regions. 

Using [61, equation 9.1.16] it can be shown by matching 
the Minkowski and de Sitter modes at t = that 



Bk 



TT UoUJkiO) 



2h 



[zFi {2t,L0k (0)) - yo(2t.Wfc (0))] ,(106) 



UoujkiO) 



2h 



'iJi{2tsUJk{Q))-Jo{2tsLOk{0))] 



(107) 



It can easily be verified using [61, equation 9.1.16] that these 
mode solutions satisfy the normalization condition (71). 



3. Particle production 

The lack of a (globally) time-like Killing vector for the de 
Sitter universe means it is not possible to unambiguously de- 
fine a Fock vacuum for all times. Additionally, the particle 
number of each mode after some expansion is observer de- 
pendent. However, for our analogue model, we circumvent 
this complication by associating an instantaneous Minkowski 
vacuum at each point in time — that is, we project into the 
quasiparticle basis that diagonalizes the many body Hamilto- 
nian to second order. This prescription we now follow has 
been outlined in Sec. IV D. Figure 1 shows the scaling factor 
for the relevant temporal regions. 

For convenience we define 



Rk{t) = -Uk{0)v ^ ^tsUkiO)e 



'tilts 



(108) 



as well as 



Ao = V2(:^(0)no/fi^fe(0). (109) 
Using the mode solutions (104) and (105) we then have 



^Irif Wl' = -L (^^-^^-W)' | \^^{Rkm + Y^{Rkm\ JoiRkit)) + [jURkiO)) + J^iRkiO))] Yo\Rk{t)) 

- 2 [Ji(i?fe(0))yi(i?fe(0)) + Jo{Rk{0))Yo(Rk{0))] Jo(i?fe(t))ro(i?fe (110) 

and 

nolxf (t)l' = A2l!l!f|#ll!| [Y,\Rkm + Y.'iRkiO))] Jl{Rk{t)) + [J?(i?fe(0)) + Jo'(^fc(0))] Y^{Rk{t)) 

- 2 [Ji(i?fc(0))yi(i?fe(0)) + Jo(i?fc(0))ro(i?fe(0))] Ji{Rk{t))Yi{Rk{t))\ . (Ill) 



We can therefore use (88) to calculate the particle production 
explicitly - for example. Fig. 2 shows the particle production 
for a fixed expansion X = 2000 with several different rates of 
expansion tg. 

4. Limits 

The particle production from the de Sitter expansion in- 
terpolates between two opposite limits; a sudden expansion 
for ts 0, and an analytically tractable asymptotic limit for 

l<2t,a'fc(0)<Xi/2. 



Sudden expansion: In this case we have the 
Minkowski mode functions: 

no\xf(.t)\'^\\l, (113) 

and the particle production from (88) reduces to the ex- 
pression (94) we found previously for a sudden expan- 
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sudden 

ts = lx 10-5 

2 X 10-5 
5 X 10-5 
1 X 10--* 
1 X 10-3 



50 100 150 200 
\k\L 



FIG. 2: Particle production for a de Sitter spacetime in the FRW 
analogue model with X — 2000 and a range of expansion rates ia as 
shown. The dimensionless scaling unit is ts — ts h/(mL'^) and the 
nonlinearity is UoNq = 10^ h'^ /{mL'^). 



sion 



1 / , , N 2 



(114) 



This behaviour can be clearly seen in Fig. 2 where the 
particle number approaches the sudden result for faster 



expansion rates (i.e., smaller values of tg). 

Asymptotic expansion: Consider the condition 

1< 2i,tjfc(0) < ^1/2, (115) 
which places a constraint on the input frequencies 

i?<a'fe(0)<i?Xi/2. (116) 

In view of the fact that ujk{t) — w/c(0)/X^/^ we can 
deduce 



^ 2UuJk{t) < 1, 



(117) 



and so rewrite the initial constraint (115) and its conse- 
quence (1 17) as 



1 « Rk{0) « 



(118) 
(119) 



In terms of the expansion X we can rewrite (88) as 



X4,na' " 
Then using (110) and (1 1 1) we have 



(120) 



7r\2 (2t,wfc(0))2 



2 J 4\/T 

- 2 [Ji(i?fe(0))ri(i?fe(0)) + Jo(i?fe(0))yo(i?fe(0))] [MRk{t))Yo{Rkit)) + MRkitWiiRkim I - i (121) 



Now using the limits for Bessel functions given in 
Abramowitz and Stegun [61], as well as some basic 
trigonometric identities, we deduce 



Nk 



(122) 



Further noting that the low k (and/or high temperature) 
limit of the Bose-Einstein distribution is 



Nk - kBT/hujk 
leads to the temperature 



(123) 



(124) 



in units of ?i = c(0) = 1. That is, the particle pro- 
duction approaches the usual cosmological result of a 
thermal spectrum [1] for an sufficiently large expansion 
that proceeds sufficiently slowly. 



VI. QUASIPARTICLE PRODUCTION BEYOND THE 
ACOUSTIC APPROXIMATION 



In this section we consider quasiparticle production where 
the acoustic approximation is not enforced — we first intro- 
duced this more general case in section lllB. In general the 
usual calculations become unmanageable in this regime, and 
so, we present only two expansion scenarios here: We first 
calculate quasiparticle production for a sudden transition in 
Sec. VIA, for which an exact solution can be easily found. In 
Sec. VI B we additionally consider the cyclic universe model. 
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for which we are able to make some general statements with 
regards to quasiparticle production. 



A. Sudden Expansion 

Particle production for a sudden expansion can be calcu- 
lated using the Bogoliubov theory from Sec. IV. The scaling 
function is given by (89) with the substitution for laboratory 
time i ^ t. The excitation of each mode is attributed to the 
quantum depletion corresponding to the initial vacuum state 
with nonlinearity Uq, projected into the Bogoliubov basis cor- 
responding to the final nonlinearity Uq/X. We can thus cal- 
culate the number of Bogoliubov quasiparticles directly using 
(57), assuming the initial state is the Bogoliubov vacuum de- 
fined by 6j^'|0) = 0. Using (88) the result is 



^out 



out in out„,in\2 



k <y 



(1 



^out2) 



(125) 



tion of the quasiparticle modes, and in particular to paramet- 
ric resonance [63]. Moreover, by taking b{tf) = b{t = 0), 
a cyclic universe is an interesting counter-point to the case of 
a sudden transition: in a sudden transition the field does not 
evolve and any particle production is entirely attributed to a 
sudden change in the effective spacetime (i.e., we project into 
a new quasiparticle basis that depends on the final nonlinearity 
U ~ Uq/X); however, in a cyclic universe model, the initial 
and final effective spacetimes are the same (so long as conden- 
sate does not evolve too far from the initial ground state) and 
particle production occurs due to parametric excitation only. 
A suitable function can be expressed as 



/2TTmt\ 



(127) 



with 6(0) = 1 and where we have defined m as the number of 
cycles of the oscillation (we note that ts = tf /mm the period 
of each oscillation), and X is now defined as the amplitude of 
the oscillation (rather than the final expansion). 



where ^Jf and is given by (60) with fj = ?7o and [/ = 



C/q/-''^ respectively. The particle production from (125) is sup- 
pressed for modes of large momenta, as can be seen from (60) 
by the observation that ^ as |k| — > oo; — this is 

consistent with the preceding discussion of Sec. IIIB. 

This calculation can also be repeated using the usual meth- 
ods of quantum field theory in curved spacetime as outlined in 
Sec. Ill — similarly to the previous calculation for the sudden 
transition in the acoustic approximation outlined in Sec. V A 
— but instead by using the field equation (52), which includes 
nonlinear dispersion. The result is [17] 



^out 



l/3k 



(126) 



This takes the same form as the result from the acoustic ap- 
proximation given by (93). However, here lo^^ and w^^* are 
the Bogoliubov mode frequencies given by (62), which in- 
cludes nonlinear dispersion, for the in and out regions respec- 
tively. 

It is straightforward to show that (125) and (126) agree ex- 
actly for all |k|. We further note that for small |k|, the Bogoli- 
ubov coefficients are Uk ~ Vk and the Bogoliubov mode ex- 
pansion coincides with the mode expansion for the quantised 
phase fluctuations (36). Therefore in this limit the acoustic 
result (94) agrees with (125). 



1. Parametric resonance 

The phenomena of parametric resonance has been previ- 
ously investigated in Bose condensed systems [64, 65], and in 
the context of an expanding universe model [66, 67]. The con- 
dition for parametric resonance occurs close to cj(fc) — Vl/2 
where Vl = 2TTm/tf is the driving frequency of some exter- 
nal parameter [63], in this case the nonlinearity. Using the 
Bogoliubov excitation spectrum e(fc) = hjj{k) from (62) then 
gives a simple estimate of the peak wave-vector for resonance: 



r, \ 1/2 



,Nof + {hTl/ts) 



2ll/2_ 



1/2 



(128) 



where t/ave = |(1 + \/X)Uq is the average nonlinearity dur- 
ing the evolution. However, it is worth noting that this result 
requires that A = 1 — 1/X from (127) is a perturbative pa- 
rameter with A <C 1; we therefore henceforth refer to this as 
the perturbative resonance condition. 

The extension of the analysis to a strongly driven system 
(A < 1) has been considered in [64]. The general result is 
that the mode frequency region where parametric resonance 
occurs broadens as A increases; that is, for the parametric 
resonance peak broadens in momentum space with larger X. 



B. Cyclic universe (2 + 1 dimensions) 

The final scenario we consider is that of a cyclic (or oscil- 
latory) universe. While models of a cyclic universe certainly 
exist in the literature [62], they are not as firmly established 
as inflationary models (such as the de Sitter universe). From a 
condensed matter point of view, an analogue model of a cyclic 
universe is interesting because it leads to parametric excita- 



VII. THE CLASSICAL FIELD METHOD 

Classical field methods are a powerful tool for approxi- 
mating the dynamics of quantum systems. Their applica- 
tion to BECs include the closely related methodologies of 
the finite temperature Gross-Pitaevskii equation [68, 69], the 
truncated Wigner approximation [11-16] and the positive-? 
method [11, 70]. 
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Neglecting quantum and thermal fluctuations, the conden- 
sate dynamics are determined by the GPE (7). This approx- 
imation arises from the assumption that the condensate is 
highly occupied so that Nq = J c?x(-0t(x)?/'(x)) ^ 1. In 
the classical field approximation this description is extended 
by also including the non-condensate modes from a low en- 
ergy subspace of the system; these modes are to be considered 
classical in that they are highly populated — this is akin to the 
Bogoliubov approximation where the commutators can be ne- 
glected. We thus proceed by expanding the field operator in 
some basis 

^(x,f) =^(/.k(x)ak(i), (129) 

k 

and replacing this in (5) by the classical field 

V'(x,i) = ^0k(x)ak(<), (130) 
kec 

where Ok ak for those modes where A'k = (akOk) ^ 1 
is satisfied; we denote these modes in the low-energy sub- 
space by fc e C. For a homogeneous Bose gas in a box with 
periodic boundary conditions the modes of the system (the 
eigenstates that diagonalize the single-particle Hamiltonian) 
are plane wave states 

</>k(x) = ^T^"""- (131^ 

The low-energy subspace C is then determined by a momen- 
tum cut-off, below which all modes are retained in the classi- 
cal field. This is formalised by the use of a projector which is 
defined by its action on some function /(x) as 

P{/(x)} = M^) f d^' 0k(x')/(x')- (132) 
kec •' 

A. Projected Gross-Pitaevskii equation 

Neglecting all modes orthogonal to C the projected Gross- 
Pitaevskii equation (PGPE) is given by 

^^d^ = V;(x) + C/oP{|^(x)p0(x)}. (133) 

For consistency with our FRW analogue model we have not 
included an external potential here. The projector is required 
for the following reasons: 

(i) The classical field approximation naturally divides the 
system into a coherent region, which is described by 
the propagation of a classical field, and an incoher- 
ent region which is neglected in the present formalism. 
In equilibrium, the system is then described by a mi- 
crocanonical ensemble since particle numbers are con- 
served by the Hamiltonian (1). 

(ii) While the finite size of the spatial grid inherently de- 



fines a momentum cut-off, the split operator Fourier 
methods used to propagate the classical field can intro- 
duce aliasing if a projector is not applied explicitly. 

(iii) The application of a projector is consistent with using a 
contact potential to describe the two-body interactions; 
this description leads to ultraviolet divergences at large 
momenta and so a cut-off is required. 

B. The truncated Wigner approximation 

A formal framework for the ideas outlined above is pro- 
vided by the truncated Wigner approximation (TWA). We 
briefly outline the method but the reader is referred to [11- 
16] for further details. 

The TWA is a phase space method originating from the 
representation of the density operator in terms of the Wigner 
function, which is familiar from quantum optics [71]. The 
master equation for the multimode density operator can be 
formally mapped to a third-order differential equation for the 
Wigner function by the application of operator correspon- 
dences. The approximation involved in the TWA is to ne- 
glect the third-order derivative terms, which become small for 
highly occupied modes. The resulting Fokker-Planck type 
equation has no diffusion term and is equivalent to evolving 
a classical field of the form (130) with the GPE, however with 
two crucial modifications: 

1. Quantum vacuum fluctuations are included in the ini- 
tial state by adding classical noise sampled from the 
Wigner distribution; the form of this noise depends on 
the Wigner function for the ground state of the system. 
For a EEC at T = the initial amplitude of each mode 
is a random Gaussian variable that is distributed accord- 
ing to the Wigner function for a coherent state. 

2. The moments of the the Wigner function give the expec- 
tation values for symmetrically ordered operators. In 
practice calculating the expectation value of an observ- 
able O requires an ensemble average over many trajec- 
tories in phase space. We denote such an expectation 
value by {0)w- 

The initial condition is given by superposition of the ground 
state and noise sampled from the Wigner distribution. 

^(x,i = 0) = V^o(x) + (5V'(x). (134) 

We can expand the noise term via a Fourier transform as 

^W^^tE^"''"'?^- (135) 

Within the truncated Wigner approximation the initial vac- 
uum state is prepared by specifying noise on each of the Bo- 
goliubov modes: 

5V(x)= (C^k(x)/3k + Vk*(x)/3£), (136) 
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where f/k(x) and Vk(x) are the plane wave modes with ampli- 
tudes Uk and respectively (as defined in Sec. IV A). The /3k 
are complex random variables that obey the Gaussian statistics 
[11,72]: 

(/3p/3q) - (/3;/3*)-0, (137) 
(/3;/3q) - (138) 

The initial state is thus constructed by populating the Bo- 
goliubov modes with half a particle per mode according to 
the TWA prescription, for the initial nonlinearity J/q; for our 
cosmological model this corresponds to the instantaneous vac- 
uum state (Minkowski vacuum) in laboratory time at t = 0. 

The Wigner and quantum expectation values for the popu- 
lation of the k mode in the Bogoliubov quasiparticle basis are 
related by 

{P^P^.)w = {{bib^}) = (iVk) + i (139) 

where the braces require that the symmetrised operator should 
be taken. The vacuum state therefore corresponds to half a 
particle per mode in the classical field. 

C. Validity of the TWA 



pressed as 

7^(x,t) = e-^^*/''[V;o(x) + <5V(x,t)], 
5V'(x,t) = ^—ZY,{ukP-i.{t)+Vk(im) (140) 

k 

where the time-dependent amplitudes are given from (57) by 

/3k(i) = WkakW ^ - Wka*-kW (141) 
|ao| |ao| 

and where the condensate phase factor is given by ao/ |q;o I = 
g-jpt/ft Referring to (139), the quasiparticle number in the 
TWA prescription (where the expectation value is implicitly 
assumed) is given by 

Nv{t) ^ {(3\{t) P^[t))w - \- (142) 

We use this quantity to calculate the quasiparticle mode pop- 
ulations in our simulations; at each time the mode functions 
Uk and Vk are determined from (59) using nonlinearity U{i). 
Hence this result is consistent with (88), as it requires pro- 
jection into the quasiparticle basis that instantaneously diago- 
nalizes (to second-order in quasiparticle operators) the many 
body Hamiltonian (1). 



In the present application of the TWA, only the k = con- 
densate mode is macroscopically occupied, the other modes 
being initially unpopulated {i.e., the quantum expectation 
value is (VVk) = 0). The requirement that iVk ^ 1 for 
each mode in the classical field is then violated. However a 
more detailed treatment of the validity of the TWA leads to 
the criterion that ':$> M for a system of N particles and M 
modes [14]. This criterion has been made explicit by Norrie 
et al. [16] as the requirement that the particle density exceeds 
the commutator for the restricted field operator . It has been 
shown that for a homogeneous system these two criteria coin- 
cide [73]. Therefore the TWA can still be applied when most 
of the modes are unoccupied as long as the average particle 
density is sufficiently large. 

Moreover the above choice for the initial state can lead to 
heating as is evident by a transient thermalization of the sys- 
tem [11, 12]. In this case the classical field dynamics deviate 
from the Bogoliubov theory valid for a weakly interacting gas; 
the system evolves to thermal equilibrium via the nonlinear 
interactions between Bogoliubov modes. This effect can be 
suppressed by evolving the classical field only for short times 
and by choosing a regime where system is weakly interacting 
{i.e., U small). This is an important consideration for our sim- 
ulations where any thermalization could obscure the effect of 
particle production. 



E. Numerical details 



Our "universe" is specified by a box with dimensions = 
jxL, Ly ~ JyL, and ~ jzL. In what follows we assume 
7a; = 7i/ = 1 and that 7^ is strictly less than one as required 
by the quasi-two-dimensional model. To facilitate numerical 
computation we introduce the dimensionless parameters 



X 



t = t 



(143) 



With a time-dependent nonlinear interaction, the two- 
dimensional GPE (35) then takes the dimensionless form 



■CNL{t)\^[' 



(144) 



We have taken Fext = for the homogeneous system. The 
effective nonlinearity is (integrating over the z direction for 
the quasi-two-dimensional geometry) 



ft2 



D. Quasiparticle number 

Following the discussion on Bogoliubov theory in Sec. IV, 
the classical field for the homogeneous system can be ex- 



Note that the wave function is normalized to unity here. The 
corresponding dimensionless speed of sound is 



m L 



NL- 



(146) 
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For completeness, we note the Bogoliubov excitation spec- 
trum (62) in dimensionless units is given by 



(147) 



In dimensionless units the spatial coordinates span the re- 
gion — i < X < i . For convenience we henceforth drop the 
bar notation (unless otherwise specified). 

Equation (144) is propagated using the 4th order Runge- 
Kutta algorithm in the interaction picture [74]. For the results 
presented here the time step was chosen so that the total nor- 
malization change during each trajectory was Anorm < 10~^ 
(for our choice of total particle number, this corresponds to a 
total loss or gain of much less than one particle for the entire 
field). 

The field is discretized on a grid of 128 x 128 points. The 
projector retains all modes with |fc| < 32 x 27r, an area which 
includes AI = 3209 modes in the classical field. 

In practice, the mode populations Afk(fc.T, ky) were resam- 
pled on polar coordinates |fc| and and then averaged over 
angle. 



F. Suitable parameter regime 

It is appropriate, at this point, to determine a viable set of 
parameters for the classical field simulations. The choice of 
simulation parameters is constrained by three main factors: 

(i) The criterion for the validity of the classical field 
method (i.e., for the TWA). 

(ii) The requirement that all modes of the system are in the 
phononic regime at the start of the simulation. In this 
regime the particle production is significant. 

(iii) A set of parameters that are experimentally relevant. 

The criteria for the validity of the TWA has been discussed 
in Sec. VllC. In particular, for our simulations, we are re- 
quired to choose a condensate population with Nq ^ M = 
3209. Additionally, noting that in our simulations the classical 
field is normalised to unity, the requirement that the system is 
weakly interacting is satisfied when the nonlinearity Catl is 
small compared with the condensate population Nq. 

We wish to investigate a regime where a significant frac- 
tion of the modes are phononic (as determined by fc < 1 /Q 
so that we can compare our results with the analytic cal- 
culations in the acoustic approximation. In computational 
units the phonon to free-particle cross-over is determined by 

2 

= Cjvl(O). That is, we require a large initial nonlinear- 
ity. 

Recent experimental observations indicate that *^Rb con- 
densates have the most widely tunable interactions via a Fes- 
hbach resonance. We refer in particular to experimental re- 
sults from the JILA group [75, 76]. In their results (see [76]) 
a stable condensate of 10^ atoms was formed with a variation 
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FIG. 3: de Sitter model: density plot for renormalized wavefunction 
at beginning (t = 0) and at end (t = tf) of expansion. Parameters 

are ts = 1 x 10"^ Civz,(f = 0) = 1 x 10^ and iVo = lO'^. 



of the scattering length from zero to 4000 oq. The associated 
diluteness factor na^ ~ 10~^ indicates such a system has sig- 
nificant interactions but can still be considered to be weakly 
interacting. 

We employ the parameters from [76], but use a larger atom 
number of iVo = 10^ while considering the same (peak) num- 
ber density Nq/V w 10^^ cm^'^. We also take the maximum 
possible initial scattering length to be a = 4000 ao at t = 0. 

With these parameters in mind we calculate the dimen- 
sionless parameters required for the classical field simulation. 
From (145) and assuming ~ "/zV^^'^ we can estimate the 
(dimensionless) initial nonlinear interaction strength is 



Cnl 



Ana 



7; 



,1/1/3 



No 



1.24 X 10^ 



(148) 



The anisotropy parameter should be taken 7^ < 1 for the 
quasi-two-dimensional geometry — we do not impose a spe- 
cific value, but note that we are free to choose a value of the 
scattering length less than 4000 qq. Therefore to meet all 
the above requirements we select CNL{t = 0) = 1 x 10^ and 
A^o = 10^ for the simulation results presented in this paper 
While a stable condensate with this atom number has not yet 
been achieved experimentally, it gives a diluteness factor less 
than na^ ^ 10^^ as found in [76], and so should in principle 
be possible. 

With a large value for A^o. we have checked that the ther- 
malization that can occur in the TWA at large nonlinearities 
is suppressed. The quasiparticle production demonstrated is 
then due solely to the effects of expanding the effective space- 
time. Additionally we note that the TWA is valid for short 
times only — this allows us to explore systems undergoing 
rapid expansion and for which there is appreciable particle 
production. 



VIII. EXPANDING UNIVERSE SIMULATIONS 

We now present the numerical results of classical field sim- 
ulations based on the TWA for the expansion scenarios out- 
lined in Sec. V — namely the de Sitter expansion, tanh and 
cyclic expansion scenarios, with the sudden transition as a 
limit of an infinitely fast de Sitter or tanh expansion. The 
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(c) ts = lx 10"'' (d) ts = lx IQ-^ 



FIG. 4: de Sitter expansion: time dependence of Bogoliubov mode populations. Parameters are CNL{t = 0) = 1 x 10^, A'o = 10^ and 
X = 2 X 10^. The blue points on each curve show where each mode crosses over from phonon to free-particle behaviour (as defined by 

2 

— Cnl)- The green dashed curve shows the analytic prediction in the acoustic approximation based on the Bogoliubov theory and 
Hamiltonian diagonalization as outlined in Sec. V A. The red solid curve shows the analytic prediction (including quantum pressure) for a 
sudden transition from (125). 



results are compared to the analytic predictions in the acous- 
tic approximation, and also to the sudden transition prediction 
that includes the nonlinear dispersion of the modes (see Eq. 
125). 

In the results shown, we have calculated the quasiparti- 
cle populations for each mode as a function of time; this 
was accomplished by projecting from the single particle ba- 
sis to the Bogoliubov basis using the expression (142), and 
using the nonlinear! ty CNL(t)- Thus the basis for counting 
quasiparticles corresponds to projecting into the instantaneous 
Minkowski vacuum at each time. 



A. de Sitter universe 

For an expansion corresponding to the de Sitter universe, 
the scaling function h{t) takes the form (100). An intuitive 
picture of the effect of quasiparticle production is demon- 
strated by Fig. 3, which gives the field density at the initial 
and final times for the case of is = 1 x 10~^. The small 
scale fluctuations given by the initial quasiparticle (Bogoli- 
ubov) vacuum are amplified to a larger scale after expansion 
has occurred. 

Figure 4 shows the results for an expansion of X = 2000 



and four different rates of expansion = 1 x 10^^,5 x 
10^^, 1 X 10^^ and 1 x lO^'^. In particular, the mode popu- 
lations are shown as a function of |fc| and time. 

For comparison, the sudden transition result for X = 2000 
from (126) is shown at the final time by the red dashed curve. 
This gives the upper limit (including the effect of the quantum 
pressure) on the permissible particle production in each mode. 
The green dashed curve shows the analytic prediction at the fi- 
nal time, that is calculated in the acoustic approximation using 
Bogoliubov theory and Hamiltonian diagonalization as out- 
lined in Sec. V A. Also shown is the time tc when each mode 
crosses from phonon to particle-like behaviour due to the ex- 
pansion of the universe, as determined by fc^/2 = CNL{tc) - 
this is indicated by the blue points on each plot. 



B. tanh Expansion 

In this case the scaling function h{t) is given by the para- 
metric curve in t determined by equations (97) and (98). The 
form of b{t) tends to exhibit a rapid early change followed by a 
long tail as it approaches its final value b{tf). Due to the com- 
putational expense in running simulations for long times, we 
took the final simulation time as — is tanh^^ (0.999). The 
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FIG. 5: tanh expansion: time dependence of Bogoliubov mode populations. Parameters are CNL{t = 0) = 1 x 10^, A'^o ~ 10^ and X = 10^. 

2 

The blue points on each curve show where each mode crosses over from phonon to free-particle behaviour (as defined by — Cnl)- The 
green dashed curve shows the analytic prediction in the acoustic approximation from (99). The red solid curve shows the analytic prediction 
(including quantum pressure) for a sudden transition from (125). 
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FIG. 6: Equipartition of energy at the final time for de Sitter expan- 
sion with four different expansion rates. Parameters are CjvL(f = 
0) = 1 X 10^ iVo = 10^ and X = 2 X 10^ 



final value of the scaling function b{tf) then reached « 99.9% 
of its target value 1/X. This minor discrepancy had a negli- 
gible effect on the results shown. 

Figure 5 shows the results for an expansion of X ~ 1000 
and two different rates of expansion, tg ^ O.l and 0.5. Sim- 
ilarly to the de Sitter expansion results, the sudden result for 
X = 1000 is shown by the red solid curve and the phonon to 
particle-like crossover is indicated by the blue points. More- 
over, the green dotted curve shows the particle production for 
each mode at the final time according to the prediction for the 
acoustic approximation given by (99). 



C. Cyclic universe 

In this case, we consider the oscillating scaling function 
given by (127). In particular, we consider two subcases: a 



single cycle (m = 1), which corresponds to a single expan- 
sion and contraction for the effective space time; and multiple 
cycles (m > 1), which correspond to m expansion and con- 
tractions. The peak wave-vector fcres for parametric resonance 
from (128) is shown by the position of the red solid line in 
each case. 



1. Single cycle (m — 1) 

We further consider two specific sets of parameters: (i) 
m = 1, X = 2000 and = 1 x 10"^; and (ii) m = 1, 
X = 2000 and = 1 x lO""'. The corresponding simulation 
results are given by Figs. 7(a) and (b). In both cases there 
is a transient phase where there is dramatic quasiparticle pro- 
duction in the time-dependent Bogoliubov basis; however, the 
net quasiparticle production at the end of the cycle (t = tf)is 
small in both cases, and in particular is very close to zero for 
the faster cycle in Fig. 7(a). Note that for = 1 x 10"'"' the 
resonance condition occurs for a wave-vector larger than the 
cut-off for the projector (/ .6.^ /i^res ^ ^cut-off)- 



2. Multiple cycles (m > 1) 

The excitation of the field due to parametric resonance is 
much greater when there are multiple cycles of the scaling fac- 
tor, and we can therefore consider a smaller driving amplitude 
X. In particular, we consider two specific sets of parameters 
for this subcase: (i) m = 5, X = 100 and = 1 x 10^**; and 
(ii) m = 10, X = 2 and tg = I x 10^''. The coiTesponding 
simulation results are given by Figs. 7(c) and (d). Clearly the 
mode population is peaked near the resonance condition. 
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(c) m = 5, X = 100, is = 1 X lO"" (d) m = 10, X = 2, = 1 X lO"" 



FIG. 7: Cyclic universe: Time dependence of Bogoliubov mode populations for four different scenarios. Parameters are C]vl(? = 0) = 
1 X 10'^, A'^o = 10^ for all cases. The red solid line shows the position of the peak wave-vector for parametric resonance from the analytic 
prediction (128). In subplots (a) and (b), the spike in population at t ~ 0.5t/ is an artifact of the time-dependent quasiparticle projection. 



rX. DISCUSSION 
A. Inflationary universe models 

The results for the de Sitter and tanh cases can be inter- 
pretted as follows. The healing length increases as the expan- 
sion proceeds; a mode k that starts as phononic {k < 1/^(0)) 
will at a later time tc cross-over to a particle-like regime 
{k ~ l/^(ic)). According to the discussion of Sec. Ill we ex- 
pect particle production to be dominant before this time, the 
mode populations becoming fixed after this time. This predic- 
tion is bourne out in the results shown in Figs. 4 and 5. 

For fast expansions, the analytic prediction in the acous- 
tic approximation overestimates the quasiparticle production 
for modes that cross-over from phonon-like to particle like. 
The role of the cross-over is further illustrated in Fig. 8. Here 
we have plotted the quasiparticle numbers Nk at four dif- 
ferent times, corresponding to the results from Fig. 4(c) for 
ts = 1 X lO^""^. For each time, the blue points represent the 
calculated from the classical field simulations, the green 
dashed curve is the analytic prediction in the acoustic approx- 
imation from Bogoliubov Hamiltonian diagonalization, and 
the solid red vertical line represents the value of the cross- 
over kc- As can be seen, for modes that are phononic (left of 
kc), the quasiparticle production from the classical field sim- 
ulations closely match the analytic prediction for the acoustic 



approximation. As the expansion proceeds kc decreases and 
Nk for modes k > kc is suppressed below the analytic predic- 
tion. 

Moreover, as is clear from Figs. 4 and 5, the sudden tran- 
sition prediction given by (126) sets an upper limit on the 
particle production in each mode. A sudden transition cor- 
responds to ^ in the de Sitter expansion case, or t , 
in the tanh expansion case; in this limit the classical field 
does not evolve, but the Bogoliubov basis that diagonalizes 
the Hamiltonian changes. The results clearly indicate that the 
particle production approaches the sudden transition predic- 
tion in both scenarios when the expansion rate is largest — 
see Figs. 4(a) and 5(a). The discrepancy for small |k| in the 
tanh expansion case is due to nonlinear interactions, which 
are neglected in the free field theory of Sec. III. We elaborate 
on this point below. 



B. Thermal equilibrium and the adiabatic regime 

There are two factors which can lead to a thermal spectrum 
for the occupation numbers of quasiparticles. 

Firstly, if the expansion is adiabatic, the particle production 
leads to a thermal spectrum [1]. This is the case even when 
the underlying theory is a free field as in Sec. II. 

Secondly, as is the case with our simulations, the CFM is 
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based on the field dynamics for the Hamiltonian (1) which im- 
pHcitly includes higher-order terms not present in the approxi- 
mate Bogoliubov Hamiltonian (61); therefore the Bogoliubov 
modes are interacting — albeit weakly — and even in the ab- 
sence of damping the system will eventually approach thermal 
equilibrium due to ergodicity. 

In the CFM the temperature for a weakly interacting system 
can be estimated by assuming equipartition of energy. Follow- 
ing Davis et al. [68, see Sec. VI], this gives: 

(^k)w - T?^, (149) 

where = ek + A is the energy for each Bogoliubov mode in 
a condensate with eigenvalue A, and ji is the chemical poten- 
tial. Here {N\^)w refers to the quasiparticle population calcu- 
lated in the Bogoliubov basis. We can then write 

Defining the dimensionless temperature as T = ksT / e^N , 
where = mL'^/h'^, and rearranging (150) gives 

That is, if the system is in thermal equilibrium the function 
Tdkl) should be constant. 

In Fig. 6 we plot this function for the de Sitter expansion 
results with = 1 x 10-^ 5 x IQ-s, 1 x IQ-"' and 1 x 10"^ 
and at the final time in each case for X = 2000. For the fastest 
expansion with = 1 x 10^^, T is approximately linear, in- 
dicating the system is not in thermal equilibrium. This is ex- 
pected since the fastest expansion rate approaches the sudden 
expansion case, for which the particle production (125) is not 
thermal. For the slowest expansion with ts — 1 x 10^^, there 
is negligible particle production so that the mode populations 
are fixed at the initial value of half a particle per mode for the 
classical field. In this case it follows from (151) that T ~ 
as evident in the plot. 

In contrast, we note for the two intermediate expansion 
rates (ts = 5 x 10~^ and 1 x lO^''), T is relatively flat for small 
|k| which corresponds to the regions in Fig. 4 where the par- 
ticle production is most significant. In these cases, the slower 
expansions result in an approximately thermal spectrum for 
the phononic modes, which is consistent for adiabatic expan- 
sion in the free-field theory. For the larger |k| modes, no par- 
ticle production occurs, and the mode populations are frozen 
at the initial value of half a particle per mode in the classical 
field. If the field was further evolved at the final nonlinearity 
CNL{t = 0)/X, the system should eventually reach thermal 
equilibrium via ergodicity. This effect is evident in the tanh 
expansion results (Fig. 5) where the nonlinearity Cnl asymp- 
totically approaches a non-zero final value. In particular, the 
nonlinear mode mixing accounts for the discrepancy between 
the analytic predictions for the free-field theory and the par- 
ticle production in the low |k| modes. This effect is more 
pronounced in Fig. 5(b) where the system evolves for a much 



longer time. 



C. Cyclic universe model 

The mode spectrum that results from the implementation of 
a cyclic universe model is markedly different to the case of 
inflationary expansion (i.e., de Sitter, tanh or sudden expan- 
sions). Specifically, there is a non-zero wave-vector at which 
the mode population peaks, which is given by the condition 
for parametric resonance as discussed in Sec. VI B. The re- 
sults for a cyclic universe are given for a single cycle by Figs. 
7(a) and (b), and for multiple cycles by Figs. 7(c) and (d). 

For the subcase of a single cycle (to = 1) there are two 
observable effects: 

(i) A peak in quasiparticle number midway through the cy- 
cle, which corresponds to the usual notion of particle 
production due to expansion as in the inflationary mod- 
els. This is an artifact of projecting the quasiparticle 
number into the time-dependent Bogoliubov basis with 
nonlinearity J7 = L/q/^ ('-S ; Hamiltonian diagonaliza- 
tion). 

(ii) The net quasiparticle number at the end of the cycle 
is determined by projecting into the Bogoliubov basis 
with nonlinearity U{tf) = U{Q). That is, the effec- 
tive spacetime is the same at the start and end of the 
cycle (providing quasiparticle production does not lead 
to appreciable depletion of the condensate). Any quasi- 
particle production is then attributed to parametric exci- 
tation, which peaks for the wave-vector corresponding 
to the resonance condition (128). For the cycle with 
a longer period (<s = 1 x 10"^), shown in Fig. 7(b), 
the resonant wave-vector is within the projected mode- 
space (i.e., fcres < ^cut-off) and there is non-zero quasi- 
particle production at the end of the cycle. In contrast, 
for the cycle of shorter period [tg = I x lO^"'), shown 
in Fig. 7(a), the resonant condition occurs for modes 
outside the projected mode-space (i.e., fci-es > ^cut-off) 
and there is negligible quasiparticle production for the 
system modes. 

Quasiparticle production is more dramatic for the case of 
multiple cycles as we can see from the results in Figs. 7(c) 
and (d). In particular, referring to the discussion in Sec. VI B, 
for the first case with m — 5 and X ~ 100, the parametric 
resonance leads to a broad peak which can be attributed to 
the large value of X. By contrast, for the second case with 
TO = 10 and X = 2 the peak is narrower due to the smaller 
value of X. A thermal component is also evident in both 
cases for low momenta modes, due to interactions between 
the quasiparticle modes. This component is expected to grow 
with longer evolution times as the modes continue to interact 
(noting that Cnl is signficant throughout the evolution). 
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FIG. 8: de Sitter expansion: Bogoliubov mode populations at four different times. Parameters are is = 1 x 10 CNhit = 0) = 1 x 10^, 
A^o = 10^ and X = 2 x 10"^ at tf. The blue points are the simulation points; the vertical solid red line shows the crossover from phonon to 

— 2 

free-particle behaviour (as defined by — Cnl)\ and the green dashed curve shows the analytic prediction in the acoustic approximation 
based on the Bogoliubov theory and Hamiltonian diagonalization as outlined in Sec. V A. 



X. CONCLUSIONS 

In summary, we have run classical field simulations for ana- 
logue models of inflationary cosmology (de Sitter and tanh 
expansions) and also for a cyclic universe model. For the in- 
flationary models the calculation of quasiparticle production 
Nk shows the following trends: 

(i) Nk is enhanced for faster expansions (small ts) and 
larger expansions (large X). In the limit of a very fast 
expansion, the results approach the sudden result from 
(125), which is expected since the field does not evolve 
for a sufficiently fast change in the nonlinearity. Quasi- 
particle production is always suppressed below the sud- 
den prediction. 

(ii) Nk is larger for small momenta, since each mode is 
strongly coupled to the effective-spacetime in this case. 
Alternatively stated, the field equation (17) becomes 
adiabatic for large momenta, so that these modes are 
not strongly excited by the expansion. Hence the classi- 
cal field simulations demonstrate the effects of Lorentz 
violation for the effective spacetime. As expected, the 
analytic predictions within the acoustic approximation 



agree favourably with the quasiparticle production for 
small momenta. 

Finally, for the case of a cyclic universe, quasiparticle pro- 
duction is entirely attributed to parametric excitation of the 
modes as there is no net change to the effective spacetime. In 
particular, parametric resonance is observed for a mode sat- 
isfying Lo{k) = Q,/2 where Q, is the driving frequency. The 
parametric resonance peak broadens in momentum space for 
larger driving amplitude A as expected from the analysis in 
[64]. 

Our calculations (both analytical and numerical) clearly in- 
dicate that quasiparticle production should occur in a number 
of different scenarios, but it is worth commenting on relevant 
experimental studies that have appeared in the literature, and 
to then suggest possible extensions to our simple model. 

A. Possible experimental implementation 

While a specific experiment corresponding to the analogue 
FRW model we have described has not yet been implemented, 
there has been significant progress in experiments which 
might ultimately lead to this goal. In particular, the two main 
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features required by our model are: (i) a homogeneous BEC 
in a box trap; and (ii) a time-varying scattering length. 

In particular, some progress has been made in experimental 
implementations of a box trap for BECs, including a square 
well potential with high barriers on an atom chip [77] and a 
novel optical trap [78]. It should be noted, however, in both 
these experiments the resulting hard-wall potential confined 
the condensate in only one dimension. On the other hand, the 
implementation of a time-varying scattering length is easily 
achievable in a number of different atomic species by making 
use of a Feshbach resonance [36, 37]. As noted in Sec. VII F, 
a promising candidate in this regard is ^^Rb which allows 
widely tunable interactions [75, 76]. However, it should be 
emphasized that there are additional complications near a Fes- 
hbach resonance from either the inelastic processes of three- 
body recombination (TBR) [79], molecule formation [80] or 
the existence of Efimov states [81]. It should be possible to 
avoid these issues in experiments by not tuning the interac- 
tions too closely to the Feshbach resonance, or by expanding 
on a time scale too fast for TBR to have an appreciable effect. 

Irrespective of the exact details for the experimental reali- 
sation of the FRW analogue model — i.e., for a homogeneous 
condensate in a box trap — we should apply either (or both) 
of the conditions: (i) that the potential at the edge of the box 
satisfies Vbox ^ C/n, or (ii) that the time-scale of expansion 
should be very small compared with the trapping frequency so 
that the condensate remains in the ground state of the trap. 

Finally, it is worth commenting that there is already exper- 
imental evidence for excitation of a condensate {i.e., quasi- 
particle production) due to a time-varying scattering length. 
In one such experiment by Claussen et al. [82], Bose con- 
densed ^''Rb atoms were subjected to an increase in scattering 
length, followed by a hold time and then a reduction in scat- 
tering length. In that work the resulting depletion of the con- 
densate increased with decreasing rise time, except for small 
hold times (ihoid < 15^s) and small rise times (trfse ^ 20^s). 
The time scale over which the scattering length was modified 
was too small for the condensate shape to adjust dynamically. 
This experiment would therefore correspond to a contraction 
and expansion of an effective spacetime, and the dependence 
of the condensate particle loss on the rise time is consistent 
with our predictions of quasiparticle production. This effect 
has been reproduced in numerical simulations by using a gen- 
eralised Gross-Pitaevskii equation [83], although it was nec- 
essary to include the effects of TBR there because of the rela- 
tively long hold times. 

B. Outlook 

The FRW analogue model we have considered is based on 
several simplifying assumptions, the two most significant be- 
ing homogeneity of the condensate and the two-dimensional 
box geometry. The preceding discussion therefore motivates 
several directions in which the formalism could be extended 
to deal with any realistic experiments that would implement a 
FRW analogue model of an expanding unverse: 



1. In an experimental implementation of the FRW ana- 
logue model, the actual trapping potential may differ 
from a two-dimensional box trap — specifically, an 
implementation will likely require a three-dimensional 
system with a non-zero potential. The extension of 
the classical field simulations from two to three dimen- 
sions is straightforward since the PGPE (133) takes the 
same form in either case. Moreover, it is straightfor- 
ward to include a realistic trap into the classical field 
simulations by specifying the non-zero potential in the 
PGPE. However, due to the significant increase in size 
of the mode space for three dimensions, these simula- 
tions would necessarily require a considerable compu- 
tational effort {i.e., running trajectories in parallel on a 
cluster of workstations). 

2. As an alternative to the condensate in a box scenario, 
we might consider the case of a BEC in a harmonic 
trap where the trapping frequency and scattering length 
are simultaneously modified in such a way that the con- 
densate density is approximately constant at the center 
of the trap. Thus the FRW analogue model we have 
considered could be approximately reproduced near the 
center of the trap. 

3. It may be necessary to include the effects of TBR to ac- 
curately describe the dynamics close to a Feshbach res- 
onance. The inclusion of TBR into the TWA has been 
previously described by Norrie et al. [84]. 

4. Finally, the presence of a thermal cloud will certainly 
affect the results of any experiment, possibly obscur- 
ing the signal of quasiparticle production due to ex- 
pansion. Finite temperature effects may be included 
by using the classical field method whereby the phase 
space is separated into a coherent region (highly occu- 
pied modes) and an incoherent region (weakly occupied 
modes). This has been formalised in terms of either the 
finite temperature GPE [68] or the stochastic GPE [72]. 

Whether or not any of the above modifications are incor- 
porated into the model, however, one should remain careful 
that the analogy to a FRW-type universe is preserved in some 
regime of interest. 



Acknowledgments 

This work was supported by Victoria University of Welling- 
ton through two PhD completion scholarships (PJ and SW). 
PJ was also supported by the New Zealand Tertiary Education 
Commission scholarship TAD- 1054. SW was also supported 
by a Victoria University Small Research Grant. CG and MV 
were supported by two Marsden Fund grants administered by 
the Royal Society of New Zealand. PJ wishes to thank Ash- 
ton Bradley, Craig Savage, Pierfrancesco Buonsante and Blair 
Blakie for useful discussions. 



25 



[1] N.D. Birrell and P.C.W. Davies, Quantum fields in curx'ed space 

(Cambridge University Press, 1982). 
[2] S.A. Fulling, Aspects of Quantum Field Theory in Curved 

Space-Time (Cambridge University Press, 1989). 
[3] S. E. Ch. Weinfurtner, arXiv:gr-qc/0404063 (2004). 
[4] P O. Fedichev and U. R. Fischer, Phys. Rev. Lett. 91, 240407 

(pages 4) (2003). 
[5] P O. Fedichev and U. R. Fischer, Phys. Rev. D 69, 064021 

(pages 9) (2004). 
[6] R O. Fedichev and U. R. Fischer, Phys. Rev. A 69, 033602 

(pages 8) (2004). 
[7] M. Uhlmann and Yan Xu and R. Schutzhold, New Journal of 

Physics 7 (2005). 
[8] C. Barcelo, S. Liberati, and M. Visser, Phys. Rev. A 68, 053613 

(pages 14) (2003). 
[9] R.A. Duine and H.T.C. Stoof, arXiv:cond-mat/0204529 (2002). 
[10] U. R. Fischer, arXiv:cond-mat/0406086 (2004). 
[11] M. Steel, M. Olsen, L. Plimak, P Drummond, S. Tan, M. Col- 

lett, D. Walls, and R. Graham, Phys. Rev. A 58, 4824 (1998). 
[12] A. Sinatra and Y. Castin and C. Lobo, Journal of Modern Optics 

47, 2629 (2000). 

[13] A. Sinatra, C. Lobo, and Y. Castin, Phys. Rev. Lett. 87, 210404 
(2001). 

[14] A. Sinatra, C. Lobo, and Y. Castin, J. Phys. B: At. Mol. Opt. 

Phys. 35, 3599 (2002). 
[15] A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, Phys. Rev. Lett. 

94, 040401 (pages 4) (2005). 
[16] A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, Phys. Rev. A 

73, 043617 (pages 18) (2006). 
[17] S. Weinfurtner, P. Jain, M. Visser, and C. W. Gardiner, Cos- 

mological particle production in rainbow spacetimes. Work in 

progress (2007). 

[18] L. J. Garay, J. R. Anglin, J. I. Cirac, and P ZoUer, Phys. Rev. 

Lett. 85, 4643 (2000). 
[19] L. J. Garay and J.R. Anglin and J.I. Cirac and P. ZoUer, Phys. 

Rev. A 63, 023611 (2001). 
[20] C. Barcelo and S. Liberati and M. Visser, Class. Quantum Gray. 

18, 1137 (2001). 

[21] M. Visser, C. Barcelo, and S. Liberati, General Relativity and 

Gravitation 34, 1719 (2002). 
[22] C. Barcelo, S. Liberati, and M. Visser, International Journal of 

Modern Physics A 18, 3735 (2003). 
[23] C. Barcelo, S. Liberati, and M. Visser, Living Rev. Rel. 8 

(2006). 

[24] M. Visser, C. Barcelo, and S. Liberati (2001), hep-th/0109033. 

[25] S. Weinfurtner, S. Liberati, and M. Visser, J. Phys. Conf. Ser. 
33, 373 (2006), gr-qc/05 12127. 

[26] W. G. Unruh, Phys. Rev. Lett. 46, 1351 (1981). 

[27] M. Visser, Acoustic propagation in fluids: An unexpected ex- 
ample of horizons?, gr-qc/931 1028 (1993). 

[28] J. F Donoghue, Phys. Rev. D 50, 3874 (1994). 

[29] S. Liberati, M. Visser, and S. Weinfurtner, Classical and Quan- 
tum Gravity 23, 3129 (2006). 

[30] S. Weinfurtner, S. Liberati, and M. Visser, J. Phys. A39, 6807 
(2006), gr-qc/05 11 105. 

[31] S. Weinfurtner, S. Liberati, and M. Visser (2006), gr- 
qc/0605121. 

[32] S. Liberati, M. Visser, and S. Weinfurtner, Phys. Rev. Lett. 96, 

151301 (2006), gr-qc/05 12139. 
[33] M. Visser and S. Weinfurtner, Massive phonon modes from 

a bec-based analog model, cond-mat/0409639 (2004), cond- 



mat/0409639. 

[34] M. Visser and S. Weinfurtner, Phys. Rev. D72, 044020 (2005), 
gr-qc/0506029. 

[35] R. Schutzhold, G. Plunien, and G. Soff, Phys. Rev. Lett. 88, 
061101 (2002). 

[36] J. M. Vogels, C. C. Tsai, R. S. Freeland, S. J. J. M. F Kokkel- 
mans, B. J. Verhaar, and D. J. Heinzen, Phys. Rev. A 56, 1067 

(1997) . 

[37] S. Inouye, M. Andrews, J. Stenger, H.-J. Miesner, D. Stamper- 
Kurn, and W. Ketterle, Nature 392, 151 (1998). 

[38] G. EUis and T. Rothman, Am. J. Phys. 61, 883 (1993). 

[39] L. P Grishchuk and Y. V. Sidorov, Phys. Rev. D 42, 3413 
(1990). 

[40] S. Liberati, M. Visser, F. Belgiorno, and D. Sciama, Physical 

Review D 61, 085023 (2000). 
[41] T. Jacobson, arXiv:gr-qc/0308048 (2003). 
[42] S. Liberati, M. Visser, F. Belgiorno, and D. W. Sciama, Phys. 

Rev. Lett. 83, 678 (1999), quant-ph/9805023. 
[43] F. Belgiorno, S. Liberati, M. Visser, and D. W. Sciama, Phys. 

Lett. A271, 308 (2000), quant-ph/9904018. 
[44] S. Winitzki, Phys. Rev. D 72, 104011 (pages 14) (2005). 
[45] W. G. Unruh and R. Schutzhold, Phys. Rev. D 68, 024008 

(pages 14) (2003). 
[46] S. Weinfurtner, A. White, and M. Visser (2007), gr-qc/07031 17. 
[47] U. R. Fischer and R. Schutzhold, Phys. Rev. A 70, 063615 

(pages 8) (2004). 
[48] D. Mattingly, Living Rev. Rel. 8, 5 (2005). 
[49] J. Martin and R. Brandenberger, Phys. Rev. D 63, 123501 

(2001). 

[50] A. Kempf and L. Lorenz (2006), gr-qc/0609123. 

[51] L. Mersini-Houghton, M. Bastero-Gil, and P. Kanti, Phys. Rev. 
D64, 043508 (2001), hep-ph/0101210. 

[52] M. Bastero-Gil and L. Mersini-Houghton, Phys. Rev. D67, 
103519 (2003), hep-th/0205271. 

[53] N. N. Bogoliubov, J. Phys. (USSR) 11, 23 (1947). 

[54] A.L. Fetter and J.D. Walecka, Quantum Theory of Many- 
Particle Systems (Dover Publications Inc., 1971, 1999, 2003). 

[55] A.A. Abrikosov and L.P. Gorkov and I.E. Dzyaloshinski, 
Methods of Quantum Field Theory in Statistical Mechanics 
(Prentice-Hall Inc., 1963). 

[56] M. Visser, Phys. Rev. A59, 427 (1999), quant-ph/9901030. 

[57] C. Bernard and A. Duncan, Annals of Physics 107, 201 (1977). 

[58] A. H. Guth, Phys. Rev. D 23, 347 (1981). 

[59] G. Gibbons and S. Hawking, Phys. Rev. D 15, 2738 (1977). 

[60] A. Lapedes, J. Math. Phys. 19, 2289 (1978). 

[61] M. Abramowitz and I. A. Stegun, Handbook of Mathemati- 
cal Functions with Formulas, Graphs, and Mathematical Ta- 
bles (Dover, New York, 1964), ninth Dover printing, tenth GPO 
printing ed., ISBN 0-486-61272-4. 

[62] P Steinhai-dt and N. Turok, Phys. Rev. D 65, 126003 (2002). 

[63] E. M. Lifshitz and L. D. Landau, Mechanics (Elsevier, 1976), 
3rd ed. 

[64] J. Garcia-RipoU, V. Perez-Garcia, and P. Torres, Phys. Rev. Lett. 
83, 1715 (1999). 

[65] M. Modugno, C. Tozzo, and F. Dalfovo, arXiv.orgxond- 

mat/0605183 (2006). 
[66] I. Zlatev, G. Huey, and P Steinhardt, Phys. Rev. D 57, 2152 

(1998) . 

[67] B. A. Bassett and S. Liberati, Physical Review D 60, 049902 

(1999) . 

[68] M. Davis, S. Morgan, and K. Burnett, Phys. Rev. A 66, 053618 



26 



(2002). 

[69] P. B. Blakie and M. J. Davis, Phys. Rev. A 72, 063608 

(pages 12) (2005). 
[70] P D. Drummond and J. F. Comey, Phys. Rev. A 60, R2661 

(1999) . 

[71] C. W. Gardiner and P. Zoller, Quantum Noise (Springer- Verlag, 

Berlin, 2004), 3rd ed. 
[72] C. W. Gardiner, J. R. Anglin, and T. 1. A. Fudge, J. Phys. B: At. 

Mol. Opt. Phys. 35, 1555 (2002). 
[73] A. A. Norrie, Ph.D. thesis. University of Otago (2005). 
[74] B. M. Caradoc-Davies, Ph.D. thesis. University of Otago 

(2000) . 

[75] J. L. Roberts, N. R. Claussen, J. R Burke, C. H. Greene, E. A. 

Cornell, and C. E. Wieman, Phys. Rev. Lett. 81, 5109 (1998). 
[76] S. L. Cornish, N. R. Claussen, J. L. Roberts, E. A. Cornell, and 

C. E. Wieman, Phys. Rev. Lett. 85, 1795 (2000). 
[77] W. Hansel and P. Hommelhoff and T. W. Hansch and J. Reichel, 

Nature 413, 498 (2001). 
[78] T. P. Meyrath and F. Schreck and J. L. Hanssen and C.-S. Chuu 

and M. G. Raizen, Phys. Rev. A 71, 041604 (pages 4) (2005). 
[79] Y. Zhang and L. Yin, Phys. Rev. A 72, 043607 (pages 4) (2005). 
[80] E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. 

Cornell, and C. E. Wieman, Nature 412, 295 (2002). 
[81] E. Braaten and H.-W. Hammer and M. Kusunoki, Phys. Rev. 

Lett. 90, 170402 (pages 4) (2003). 
[82] N. R. Claussen and E. A. Donley and S. T. Thompson and C. E. 



Wieman, Phys. Rev. Lett. 89, 010401 (2002). 

[83] R. A. Duine and H. T. C. Stoof, Phys. Rev. A 68, 013602 
(pages 17) (2003). 

[84] A. A. Norrie and R. J. Ballagh and C. W. Gardiner and A. S. 
Bradley, Phys. Rev. A 73, 043618 (2006). 

[85] There is an assumption here, that the linear operator U is invert- 
ible, which is true only if the null space of U is the zero vector 
space. This is certainly true for the homogeneous BEC as we 
will show in Sec. IV. 

[86] This is because the projected Gross-Pitaevskii equation — 
which the classical field method is based on — takes the same 
form for d = 2 or d = 3. While the extension of the classical 
field method to d = 3 spatial dimensions is straightforward, 
this requires a larger mode space and our simulations would 
have required an unfeasibly large computational time. 

[87] The models considered in [6, 7] correspond to the case of a 
physical expansion, which corresponds to the free expansion of 
a BEC. The scaling factor there has a different physical signifi- 
cance from the scaling factor for the effective expansion that is 
considered here. 

[88] With no loss of generality (as the Bessel equation is even in rj), 
the arguments of Ji and Yi are taken to be negative with respect 
to rj but positive overall in the region of interest - this facilitates 
computation as the Bessel functions are then real quantities. 



